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Abstract 

It  has  been  shown  theoretically  that  under  mild  conditions  multidimensional  signals  can 
be  recovered  from  one-level  crossings  (e.g.  zero  crossings).  However,  the  accuracy  with  which 
locations  of  the  one-level  crossings  need  to  be  specified  is  large  enough  to  limit  the  applicability 
of  such  a  method  in  many  practical  situations.  In  this  thesis,  we  will  propose  two  major 
sampling  strategies  for  reconstruction  of  signals  from  multiple-level  crossings. 

In  our  first  approach,  we  extend  new  theoretical  results  in  multivariate  polynomial  interpo¬ 
lation  theory,  in  order  to  define  a  variety  of  semi-implicit  sampling  strategies.  These  strategies, 
which  provide  sufficient  conditions  for  recovery  of  multidimensional  signals  from  non-uniform 
samples  on  lines  of  rational  slope,  are  ultimately  applied  to  the  problem  of  reconstruction  from 
multiple-level  crossings.  Although  these  semi-implicit  results  are  general  enough  to  be  used 
for  recovery  from  signal  crossings  with  arbitrary  functions,  they  do  not  provide  conditions  for 
reconstruction  of  signals  from  an  arbitrarily  small  number  of  thresholds.  In  order  to  circumvent 
this  difficulty,  we  propose  a  second  approach  which  is  implicit,  and  uses  algebraic  geometric  con¬ 
cepts  to  find  conditions  under  which  a  signal  is  almost  always  reconstructible  from  its  multilevel 
threshold  crossings. 

A  problem  distinct  from  that  of  uniquely  specifying  signals  with  level  crossings  is  that  of 
developing  specific  algorithms  for  recovering  them  from  level  crossing  information,  once  it  is 
known  that  the  signals  satisfy  the  appropriate  constraints.  We  propose  a  variety  of  reconstruc¬ 
tion  algorithms  for  each  of  our  two  approaches,  and  demonstrate  results  for  several  images. 
Having  proposed  a  variety  of  sampling  and  reconstruction  strategies,  we  then  present  a  prelim¬ 
inary  investigation  of  their  quantization  characteristics.  In  doing  so,  we  find  that  the  dynamic 
range  and  bandwidth  requirements  for  representation  of  signals  via  multiple  level  threshold 
crossings  lie  in  between  those  of  Nyquist  and  zero  crossing  representation.  Moreover,  under 
certain  circumstances,  our  semi-implicit  and  implicit  sampling  strategies  become  identical  to 
Nyquist  sampling.  This  will  bridge  the  gap  between  explicit,  semi-implicit,  and  implicit  sam¬ 
pling  strategies,  unify  seemingly  unrelated  sampling  schemes,  and  provide  us  with  a  spectrum 
of  sampling  schemes  for  multidimensional  signals. 
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Chapter  1 

Introduction 

Signal  reconstruction  from  partial  information  has  been  an  active  area  of  research  for  many 
years.  Previous  work  in  this  area  has  involved  developing  conditions  under  which  both  one¬ 
dimensional  and  two-dimensional  signals  are  uniquely  specified  with  Fourier  transform  magni¬ 
tude,  phase,  or  signed-magnitude  information  [3, 4, 5, 6, 7].  Recently,  Curtis  et.al.  considered  the 
problem  of  signal  reconstruction  from  Fourier  transform  sign  information  for  multidimensional 
signals  [8j.  By  exploiting  the  duality  between  space  and  frequency  domains,  they  applied  their 
result  to  the  problem  of  reconstruction  from  zero-crossing  information  [9] .  These  zero  crossing 
results  are  much  less  restrictive  and  more  broadly  applicable  than  the  earlier  ones  based  on 
two-dimensional  extensions  of  one-dimensional  results. 

Representing  signals  with  their  zero  crossings  has  important  practical  and  theoretical  im¬ 
plications.  From  a  practical  point  of  view,  the  zero-crossing  results  have  applications  in  image 
processing  and  vision,  where  the  information  contained  in  the  edges  of  objects  is  considered 
to  be  important  [10],  Also,  in  situations  where  an  image  is  clipped  or  otherwise  distorted  in 
such  a  way  as  to  preserve  zero-crossing  or  level-crossing  information,  it  is  possible,  at  least  in 
principle,  to  recover  the  signal  from  its  distorted  version. 

From  a  theoretical  point  of  view,  reconstruction  from  zero  crossings  is  an  example  of  implicit 
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sampling.  Representation  of  a  function  by  its  samples  or  derivatives  taken  at  preselected  instants 
of  time  has  received  extensive  attention  in  theory  as  well  as  in  engineering  practice.  Sampling 
methods  based  on  prespecified  sampling  times  are  referred  to  as  explicit ,  in  order  to  distinguish 
them  from  implicit  methods  of  sampling,  in  which  a  representation  is  sought  in  terms  of  the 
instants  in  which  the  function  assumes  predetermined  values.  The  idea  of  implicit  sampling  was 
introduced  by  Bond  and  Cahn  [11],  who  considered  representation  and  manipulation  of  one¬ 
dimensional  signals  by  means  of  their  real  and  complex  zeros.  Extensive  work  on  this  topic  has 
been  done  by  Voelcker  [12],  and  computer  simulation  of  these  results  has  been  reported  by  Sekey 
[13].  In  addition,  Bar-David  [14]  considered  the  important  case  of  one-dimensional  implicit 
sampling  in  terms  of  real  variables  alone.  These  results  have  been  used  to  overcome  distortions 
that  are  incurred  either  by  intentional  nonlinear  processing  or  by  inadvertent  nonlinearities. 
Such  nonlinearities  might  arise  in  single  side-band  systems,  where  the  audio  signals  are  typically 
hardlimited  in  order  to  decrease  their  dynamic  range.  Also,  in  magnetic  tape  recording,  a  strong 
1  'her  frequency  bias  tone  is  usually  added  onto  the  signal  to  ensure  fidelity  in  the  presence  of 
inherent  material  nonlinearity. 

The  majority  of  the  research  in  implicit  sampling  has  been  restricted  to  one-dimensional 
problems.  Indeed,  the  zero  crossings  results  of  Rotem  and  Zeevi  [15]  and  Curtis  et.al.  [8,9] 
are  the  only  examples  of  implicit  sampling  for  multidimensional  signals.  Rotem  and  Zeevi’s 
results  are  an  extension  of  Logan’s  [16]  one-dimensional  result,  which  only  deals  with  bandpass 
signals.  The  results  due  to  Curtis  et.al.  however,  are  truly  two-dimensional,  since  they  take 
advantage  of  the  fact  that,  in  contrast  to  the  one-dimensional  case,  the  zero  crossing  contours  in 
two  or  higher  dimensions  contain  infinitely  many  points.  These  results  are  much  less  restrictive 
and  more  general  than  those  based  on  the  extension  of  one-dimensional  results.  However,  the 


accuracy  with  which  the  locations  of  the  zero  crossings  need  to  be  specified  is  large  enough 
to  limit  their  applicability  in  many  practical  situations.  In  effect,  by  representing  the  two- 
dimensional  signal  with  zero  crossings  or  threshold  crossings,  the  amplitude  information  of  the 
original  signal  is  embedded  in  the  exact  location  of  the  threshold  crossings.  Consequently,  while 
the  original  signal  can  be  sampled  at  the  Nyquist  rate,  the  threshold  crossing  representation 
may  require  a  considerably  higher,  possibly  infinite  sampling  rate  to  preserve  the  threshold 
crossing  locations  adequately.  Thus,  the  total  number  of  bits  or  the  bandwidth  required  in  the 
threshold  crossing  representation,  is  much  higher  than  that  required  by  sampling  and  quantizing 
the  original  signal.  Therefore,  the  results  on  signal  reconstruction  from  threshold  crossings  are 
more  useful  in  applications  in  which  the  exact  threshold  crossing  points  are  available. 

It  is  possible,  however,  to  view  the  representation  of  signals  with  threshold  crossings  as  a 
trade-off  between  bandwidth  and  dynamic  range,  in  the  sense  that  if  the  available  bandwidth  is 
sufficient  to  preserve  the  threshold  crossings  accurately,  then  the  dynamic  range  requirements 
are  significantly  reduced.  On  the  other  hand,  representation  of  signals  via  their  samples  at  the 
Nyquist  rate  can  be  considered  as  requiring  minimal  bandwidth  and  infinite  dynamic  range. 
This  is  because  exact  recovery  of  signals  via  Nyquist  sampling,  requires  amplitude  information 
at  prespecified  points,  to  infinite  precision.  Thus,  the  natural  question  which  arises  is  whether 
or  not  there  are  intermediate  sampling  and  reconstruction  schemes,  whose  characteristics  lie 
between  these  two  extremes.  Our  objective  in  this  thesis  has  been  to  derive  sampling  schemes 
which  bridge  the  gap  between  Nyquist  sampling  and  zero-crossing  representations  by  enabling 
us  to  recover  signals  from  multiple  level  threshold  crossings.  To  this  end,  we  have  proposed  a 
variety  of  semi-implicit  and  implicit  sampling  strategies  in  Chapters  2  and  4  respectively.  Semi- 
implicit  samples  of  a  multidimensional  signal  are  defined  to  be  points  whose  coordinates  are 


mathematically  related  to  each  other.  As  we  will  see,  some  of  our  results  on  reconstruction  from 
level  crossings  arc  general  enough  to  be  used  for  recovery  of  signals  from  their  crossings  with 
arbitrary  functions.  Possible  applications  of  these  results  are  for  conversion  of  half  tone  images 
to  continuous  tone  ones  [17],  In  addition,  our  results  can  also  be  applied  to  the  more  general 
problem  of  reconstruction  from  non-uniformly  spaced  samples.  Reconstruction  of  functions  from 
their  samples  on  nonuniformly  distributed  locations  is  an  important  task  in  many  applications 
such  as  machine  vision  [18,19],  radio  astronomy  [20],  and  computed  tomography  [21],  as  well 
as  natural  sciences  such  as  geology,  meteorology,  and  oceanography.  Several  techniques  based 
on  non-harmonic  Fourier  series  have  been  proposed  [22,23]  for  reconstruction  of  bandlimited, 
one-dimensional  functions,  sampled  at  irregular  intervals.  These  methods,  however,  are  limited 
to  sample  sequences  which  have  only  minor  deviations  from  uniform.  Methods  that  have  been 
used  to  perform  such  reconstructions  for  multidimensional  functions  include  nearest-neighbor 
and  bilinear  interpolation,  surface  functional  minimization  by  relaxation  or  gradient  descent 
methods,  and  Gaussian  smearing  resampling.  However,  these  methods  either  do  not  result  in  a 
minimum  possible  reconstruction  error,  or  they  require  a  priori  knowledge  about  the  form  of  the 
function.  In  fact,  the  only  exact  reconstruction  strategy  for  the  multidimensional  case,  proposed 
by  Clark  [24],  is  somewhat  restrictive  and  its  corresponding  algorithm  is  rather  heuristic.  As 
we  will  see,  part  of  our  results  on  reconstruction  from  level  crossings  are  general  enough  for 
recovery  of  signals  from  their  non-uniform  samples. 

In  summary,  our  objective  in  this  thesis  has  been  to  derive  semi-implicit  and  implicit  sam¬ 
pling  schemes,  whose  characteristics  lie  in  between  the  Nyquist  and  zero-crossing  representa¬ 
tions.  These  strategies  result  in  methods  of  reconstructing  multidimensional  signals  from  their 
multiple  level  threshold  crossings,  thus,  providing  a  bridge  between  Nyquist  reconstruction  from 


explicit  samples  at  prespecified  points  and  reconstruction  via  one-level  crossings.  As  it  turns 
out,  some  of  these  sampling  techniques  are  general  enough  to  be  used  in  problems,  such  as 
reconstruction  from  crossings  with  arbitrary  functions  or  from  non-uniformly  spaced  samples. 

The  outline  of  this  thesis  is  as  follows.  We  shall  begin  by  briefly  reviewing  related  re¬ 
search  on  reconstruction  from  zero-crossings  and  formulating  the  problem  of  reconstruction 
from  multi-level  threshold  crossings.  In  Chapter  2  we  develop  new  results  in  bivariate  polyno¬ 
mial  interpolation  theory  which  are  ultimately  used  to  derive  semi-implicit  sampling  strategies 
for  bandlimited,  periodic  (BLP)  signals.  Chapter  3  describes  various  reconstruction  algorithms 
for  the  sampling  strategies  of  Chapter  2.  In  Chapter  4,  we  propose  a  second  strategy  for  re¬ 
construction  from  multiple  level  threshold  crossings.  Although  the  semi-implicit  approach  of 
Chapter  2  can  be  applied  to  more  general  problems  such  as  reconstruction  from  non-uniformly 
spaced  samples  or  reconstruction  from  arbitrary  function  crossings,  the  results  in  Chapter  4 
are  less  restrictive  in  the  sense  that  they  enable  us  to  reconstruct  signals  from  their  crossings 
with  arbitrary  number  of  thresholds.  In  Chapter  5,  we  will  present  the  results  of  a  prelim¬ 
inary  investigation  of  the  quantization  properties  of  the  various  sampling  and  reconstruction 
strategies  described  throughout  the  thesis.  More  specifically,  we  will  show  how  quantization 
characteristics  of  reconstruction  from  multiple  level  crossings  vary  as  a  function  of  the  number 
of  thresholds.  As  we  will  see,  the  representation  of  two-dimensional  signals,  via  their  amplitude 
quantized  explicit  samples,  is  intimately  related  to  their  position  quantized  implicit  or  semi- 
implicit  samples.  Indeed,  sampling  and  reconstruction  from  multiple  level  crossings,  together 
with  Nyquist  and  zero  crossing  sampling,  provide  us  with  a  wide  spectrum  of  signal  representa¬ 
tion  with  different  bandwidth  and  dynamic  range  requirements.  Finally,  conclusions  and  future 
directions  of  research  are  included  in  Chapter  6. 


In  the  remaining  part  of  this  chapter,  we  will  briefly  review  the  existing  results  on  recon¬ 
struction  from  one-level  crossings,  and  formulate  the  problem  of  reconstruction  from  multiple 
level  threshold  crossings. 

1.1  Previous  Results  and  Problem  Formulation 

There  has  been  a  great  deal  of  interest  in  the  zero  crossing  representation  of  signals  in  recent 
years.  Most  of  the  results  on  the  unique  specification  of  one-dimensionai  signals  are  based  upon 
the  fact  that  a  bandlimited  function  is  entire  and  is,  thus,  uniquely  specified  by  its  real  and 
complex  zeros  to  within  a  constant  and  exponential  factor.  An  arbitrary  bandlimited  function 
is  uniquely  specified  by  its  real  zero  crossings,  if  all  of  its  zeros  are  real.  Thus,  a  number  of 
previous  research  efforts  concentrated  on  identifying  the  conditions  under  which  signals  have 
only  real  zeros  and  on  developing  methods  for  modifying  a  signal  so  that  all  of  its  zeros  become 
real  [14].  Despite  this,  most  one-dimensional,  bandlimited  signals  encountered  in  practice,  do 
not  satisfy  the  constraints  associated  with  these  results.  They  are  not  uniquely  specified  by 
their  zero  crossings  unless  they  satisfy  some  additional  constraints,  which  effectively  guarantee 
that  they  contain  a  sufficient  number  of  zero  crossings. 

Although  a  considerable  amount  of  theoretical  work  has  been  devoted  to  the  problem  of 
reconstruction  of  one-dimensional  signals  from  zero  crossings,  much  less  work  has  been  devoted 
to  the  corresponding  two-dimensional  problem.  Logan’s  result  has  been  extended  to  two  dimen¬ 
sions  [10,15]  by  requiring  a  one-dimensional  signal  derived  from  the  original  two-dimensional 
one  to  satisfy  the  constraints  of  Logan’s  theorem.  However,  the  two-dimensional  problem  is  fun¬ 
damentally  different  from  the  one-dimensional  one,  since,  in  two  dimensions,  the  zero  crossings 
are,  in  general,  contours  rather  than  isolated  points.  Curtis  and  Oppenheim  [9]  were  first  to  take 


advantage  of  this  fact;  their  main  result  states  that  for  a  BLP  signal  with  a  {2N  +  1)  x  (2N  +  1) 
region  of  support  in  the  Fourier  domain,  16 N 2  +  1  or  more  samples  of  the  zero  crossings  are 
sufficient  to  reconstruct  the  signal  to  within  a  scale  factor. 

As  shown  in  [9],  in  carrying  out  the  above  reconstruction,  the  locations  of  the  zero  crossings 
need  to  be  specified  accurate  to  16  digits.  Less  accurate  specification  of  the  crossings  results 
in  unsuccessful  reconstruction.  Thus,  by  representing  the  two-dimensional  signal  with  zero 
crossings  or  threshold  crossings,  the  amplitude  information  in  the  original  signal  is  embedded 
in  the  exact  location  of  the  threshold  crossings.  Consequently,  while  the  original  signal  can  be 
sampled  at  the  Nyquist  rate  with  many  bits  for  amplitude  specification,  the  threshold  crossing 
representation  requires  a  considerably  higher  sampling  rate  to  preserve  the  threshold  crossing 
locations  adequately.  Thus,  it  is  possible  to  view  the  representation  of  signals  with  threshold 
crossings  as  a  trade-off  between  bandwidth  and  dynamic  range,  in  the  sense  that  if  the  available 
bandwidth  is  sufficient  to  preserve  the  threshold  crossings  accurately,  then  the  dynamic  range 
requirements  are  significantly  reduced.  On  the  other  hand,  representation  of  signals  via  their 
samples  at  the  Nyquist  rate  can  be  considered  as  one  which  requires  minimal  bandwidth  and 
large  dynamic  range.  Our  main  goal  in  this  thesis  is  to  develop  intermediate  sampling  schemes 
for  reconstruction  from  multiple  level  threshold  crossings,  so  that  their  bandwidth  and  dynamic 
range  requirements  lie  in  between  Nyquist  and  one-level  crossing  representation. 

Our  approach  to  the  above  problem  has  been  to  represent  BLP  signals  in  terms  of  polyno¬ 
mials.  The  reasons  for  doing  so  are  twofold.  First,  BLP  signals,  which  represent  a  fairly  large 
and  general  class  of  signals,  can  be  written  as  polynomials,  via  their  Fourier  series  expansions. 
Second,  since  reconstruction  from  multiple  level  crossings  is  a  special  case  of  reconstruction 
from  non-uniform  samples,  we  hope  to  be  able  to  use  a  variety  of  mathematical  results  on 


polynomial  interpolation  theory.  Consider  a  BLP  signal 


/(*,»)=  E  E  nfci,M^2,r(M  +  t3y)  (ii) 

*i  =  -Af  ki  =  -N 

The  polynomial  representation  of  the  above  signal  is  given  by: 

g{WuW2)  =  /(x,y)W*W* 

2  N  2  N  (1.2) 

=  E  E  N,k2-  N)W^W^ 

ki=  0  *j=0 

where 


W j  =  e^2” 

W2  =  e;2TW 

Since  reconstruction  of  /(*,  y)  is  equivalent  to  finding  the  coefficients  of  the  polynomial  y(Wi,  Wj), 
results  from  multivariate  polynomial  interpolation  theory  can  be  applied  directly  to  a  variety 
of  multidimensional  reconstruction  problems. 

Unlike  the  univariate  case,  interpolation  with  multivariate  polynomials  is  a  non-trivial  task. 
Whereas  n  arbitrary  samples  of  an  univariate  (one-dimensional)  polynomial  of  degree  n  -  1  are 
sufficient  to  find  its  coefficients,  the  analogous  result  in  dimensions  higher  than  one  does  not 
hold,  primarily  because  there  are  no  Chebychef  systems  in  Ra  for  s  >  2.  Chebychef  systems  are 
important  in  interpolation  theory,  and  have  been  studied  extensively  by  many  people  including 
Karlin  (26),  Karlin  and  Studen  [27],  and  Krein  [28].  The  definition  is  as  follows: 


Definition  1.1  A  linearly  independent  set  of  continuous  functions  (uo(x), ...,  u„(x)}  defined 
on  [a,  6]  is  a  Chebychef  system  if  for  any  a  <  xq  <  ...  <  xn  <  b  and  Sto,...,yn  €  R,  there  is  a 
unique  linear  combination 

n 

«(*)  =  E°ju*(*) 

;= o 


satisfying 


«(*»)  =  y i 


i  =  0,  ...,n 
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Clearly,  the  set  of  n  continuous  functions  consisting  of  the  powers  of  x  form  a  Chebychef 
system.  Another  example  of  a  Chebychef  system  is 


Ui(x)  =  e 


i  =  0, 1, ...,  n 


where  A,  are  distinct  and  x  €  (-oo,+oo)  [29j.  Chebychef  systems  are  helpful  for  studying 
univariate  interpolation.  Unfortunately,  as  soon  as  we  turn  to  multivariate  interpolation,  we 
must  leave  them  behind,  since  there  is  no  set  of  n  universal  functions  which  can  be  used  for 
interpolation  at  any  n  distinct  points.  An  implication  of  this  result  is  that  powers  of  x  or  y 
do  not  form  a  Chebychef  system  in  R2 ,  and  thus,  bivariate  polynomials  are  not,  in  general, 
uniquely  reconstructible  from  their  samples  at  arbitrary  locations.  Indeed,  if  arbitrary  samples 
of  a  bivariate  polynomial  uniquely  specified  its  coefficients,  the  problem  of  reconstruction  from 
multiple  level  threshold  crossings  would  have  been  trivial.  This  is  because  samples  of  our 
threshold  contours  could  then  be  used  to  find  the  the  coefficients  of  W2)  or  equivalently 

/(*>  y)- 

There  are  two  ways  out  of  this  dilemma.  One  would  be  to  study  conditionally  regular 
interpolation  methods.  An  interpolation  method  is  called  regular,  if  it  is  uniquely  solvable 
for  any  selection  of  the  points  of  interpolation.  Conditionally  regular  interpolation  methods 
are  not  solvable  for  all  selection  of  points,  but  only  for  most  of  them  [31].  Roughly  speaking, 
they  are  uniquely  solvable  with  probability  1.  For  methods  of  this  type,  if  one  has  a  concrete 
problem  and  selects  the  interpolation  points  at  random,  it  will  be  extremely  unlikely  that  the 
problem  will  be  unsolvable.  We  will  take  this  approach  to  find  the  sufficient  conditions  for 
reconstruction  from  multilevel  threshold  crossings  in  Chapter  4.  The  second  approach  is  to 
impose  certain  restrictions  on  the  locations  of  the  interpolation  points.  These  conditions  will 
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guarantee  that  the  resulting  interpolation  problem  has  a  unique  solution.  Our  approach  in 
Chapter  2  consists  of  developing  a  series  of  conditions  under  which  bivariate  polynomials,  or 
equivalently,  two-dimensional  signals  can  be  uniquely  recovered. 


Chapter  2 

Semi-Implicit  Sampling  Approach 
to  Reconstruction 

As  we  saw  in  secton  (1.1),  the  polynomial  representation  of  multidimensional  signals  sug¬ 
gests  that  a  wide  variety  of  multidimensional  reconstruction  problems  can  be  approached  via 
results  from  multivariate  polynomial  interpolation  theory.  In  addition,  we  saw  that  the  inherent 
difficulty  in  multivariate  interpolation  is  the  lack  of  existence  of  Chebychef  systems  in  dimen¬ 
sions  higher  than  one.  In  Chapter  4,  we  will  propose  conditionally  regular  interpolation  as  a 
possible  way  to  circumvent  this  difficulty.  In  this  chapter,  however,  our  approach  is  to  derive 
a  number  of  theoretical  results  on  multivariate  polynomial  interpolation  theory  by  imposing 
constraints  on  the  location  of  interpolation  points.  These  results  will  ultimately  be  used  to 
develop  a  spectrum  of  sampling  strategies  for  reconstruction  of  multidimensional  signals  from 
non-uniformly  spaced  samples  in  general,  and  multiple  level  crossings  in  particular. 

We  will  begin  this  chapter  with  a  brief  review  of  existing  results  in  bivariate  polynomial 
interpolation.  In  section  (2.2),  we  will  derive  three  results  in  multivariate  interpolation  theory, 
which  become  progressively  more  general.  Utilizing  polynomial  representation  of  BLP  signals, 
as  described  in  the  previous  chapter,  we  apply  our  new  theoretical  results  on  multivariate  inter¬ 
polation  to  derive  a  variety  of  semi-implicit  sampling  strategies  for  the  reconstruction  of  signals 
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from  their  non-uniformly  spaced  samples  on  lines  of  rational  slope.  As  we  will  see  in  section 
(2.3),  these  semi-implicit  sampling  schemes  can  be  applied  to  the  problem  of  reconstruction 
from  level  crossings  and  to  a  variety  of  other  problems  such  as  reconstruction  from  crossings 
with  arbitrary  functions  and  reconstruction  from  projections. 

2.1  Review  of  Bivariate  Polynomial  Interpolation  Theory 

In  this  section,  we  will  review  some  of  the  existing  results  from  bivariate  polynomial  interpo¬ 
lation  theory.  We  will  be  concerned  primarily  with  interpolation  using  only  a  function’s  values 
as  opposed  to  its  derivatives.  As  we  mentioned  earlier,  our  approach  in  this  chapter  involves 
constraining  the  locations  of  the  interpolation  points,  in  order  to  guarantee  unique  recovery  of 
the  polynomials  and  their  associated  BLP  signals. 

Bivariate  polynomial  interpolation  can  be  done  either  in  Il„,  the  space  of  polynomials  with 
total  degree  less  than  or  equal  to  n,  or  in  Il(n  m),  the  space  of  polynomials  p(w,  z)  with  maximum 
degree  n  in  w  and  m  in  z.  The  total  degree  of  a  polynomial  p(xv,z)  is  defined  to  be  the 
degree  of  the  one-dimensional  polynomial  p(w,w).  The  regions  of  support  of  the  coefficients 
of  polynomials  in  rij  and  11(2,2)  are  illustrated  in  figure  (2.1).  The  most  general  result  on 
interpolation  in  ITn  was  derived  by  Gasca  and  Maeztu  [36)  in  1982.  The  basic  idea  behind 
their  results  is  to  choose  a  set  of  straight  lines  r*  in  if2,  each  of  which  is  associated  with  a 
polynomial  of  first  degree  in  w  and  z.  With  each  line  rt,  we  consider  a  set  of  straight  lines 
in  such  a  way  that  the  intersections  determined  by  on  r,  are  the  points  at  which  the 
interpolation  data  will  be  given.  The  lines  r,  and/or  may  appear  with  multiplicity  greater 
than  one,  leading  to  derivative  values  as  interpolation  data.  In  general,  the  formulation  with 
derivatives  results  in  the  Hermite  interpolation  problem,  although  the  simplest  case  reduces  to 


Lagrange  interpolation  which  only  uses  the  sample  values.  In  any  case,  products  of  lines  r,  and 
rj'*  are  constructed  in  order  to  obtain  a  set  of  polynomials  spanning  the  vector  space  in  which 
the  interpolation  problem  has  a  unique  solution  (36).  The  most  important  consequence  of  the 
theoretical  developments  of  Gasca  and  Maetzu  can  be  stated  els  follows: 

Theorem  2.1  (Gasca  and  Maetzu  [36])  Consider  the  distinct  lines  lo,  with  the  set  of 
distinct  points 

on  /,.  If  none  of  the  interpolation  points  are  on  the  intersection  of  any  two  lines  from 

/0). then  for  any  data  set 

{tjl)  I  3  =  0, »  =  0, 1, ...,  n} 
there  is  a  unique  bivariate  polynomial  p  €  n„  such  that 

p{^\  *}*’)  =  <S$)  0  <j<i,  o  <j<n 

The  proof  is  included  in  [36,29,30],  and  an  example  of  the  geometric  distribution  of  the 
sampling  points  required  by  this  theorem  for  n  =  2  is  shown  in  figure  (2.2).  A  special  case  of 
the  above  theorem  was  proposed  earlier  by  Stenger  [35J  and  Chung  [32].  Their  results  primarily 
deal  with  the  case  where  the  interpolation  points  are  chosen  on  parallel  lines.  The  proof  for  the 
case  where  the  interpolation  lines  are  chosen  parallel  to  the  x  axis  is  rather  straightforward, 
and  can  be  briefly  outlined  in  the  following  manner.  Let  us  denote  the  points  on  the  ith  line 
parallel  to  the  x  axis  by 

{(«>,,*)  I  j  =  0,l,...,i} 

We  cam  rewrite  the  polynomial  p(w,  z)  as: 


v  ^3 


.y 


p{w,z)  =  qn(z)  +  wqn-i(z)  +  ...  +  wn  1qi{z) 


21 


(2.1) 
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Figure  2.2:  Geometric  distribution  of  the  s&mpling  points  of 
Theorem  (2.1)  for  n-2. 

where  9t(z)  is  a  polynomial  of  degree  k  —  1  in  z.  In  the  first  stage,  we  substitute  the  n  +  1  points 
on  l„  into  the  above  equation  in  order  to  find  qk(zn)  for  1  <  k  <  n.  Since  <71(2)  is  a  constant, 
its  value  can  be  determined  by  qi(zn).  In  the  second  stage,  the  n  points  on  /n_j  are  used  to 
specify  qk(zn-i)  uniquely  for  2  <  A:  <  n.  Using  the  value  of  <72(2,,)  from  the  first  stage  and 


?2(*n-i)  from  the  second  stage,  we  can  uniquely  specify  the  two  coefficients  of  the  polynomial 


92(2).  Similarly  in  the  ith  stage  we  first  use  the  points  on  the  (n  +  1  -  f)th  line  to  find  qk(z) 


for  1  <  k  <  i  —  1,  and  then  use  qi{zi)  for  n  +  1  —  i‘</<nin  order  to  find  the  1  coefficients  of 


qi{z).  Repeating  this  procedure,  we  can  uniquely  determine  the  polynomial  <7,(2)  for  1  <  1  <  n, 


so  that  p(u;,z)  becomes  completely  and  uniquely  determined. 


Most  of  the  results  in  bivariate  polynomial  interpolation  theory  deal  with  interpolation  in 
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Iln,  the  space  of  polynomials  of  total  degree  less  than  or  equal  to  n.  The  only  existing  result 
on  interpolation  in  FI(n  m),  the  space  of  polynomials  p(u;,  z)  of  maximum  degree  n  in  w  and  m 
in  z,  deals  with  the  case  where  the  interpolation  points  are  on  a  nonuniform  rectangular  grid 
and  can  be  stated  as  follows: 


Theorem  2.2  (Non-uniform  Rectangular  Sampling  [29,32,33,34])  Given  a  set  of  points 

{(«;,,  z>)|t  =  0,  j  =  0,  } 

and  data 

{ttj\i  =  0,...,n;  j  =  0,  } 

there  exists  a  unique  bivariate  polynomial 


p(w,z)  =  5^X]a(«,i)u;V 
i=0j=0 


such  that 


P(w«.2j)  =  Uj 


i  =  0,  n;  j  =  0, m; 


lV 


The  proof  is  straightforward  and  is  given  in  many  references  including  [29,32,33,34].  In¬ 
tuitively,  for  a  fixed  value  of  j  =  jo,  the  points  (u/0,  zy„), (u>n,  zju)  can  be  used  to  find  the 
coefficients  a(0,  jo), ...,  a(n,  jq).  Repeating  this  procedure,  we  can  find  all  the  coefficients  of 
p(w,z).  An  example  of  the  geometric  distribution  of  the  sampling  points  required  by  this  the¬ 
orem  for  n  =  2  and  m  =  1  is  shown  in  figure  (2.3).  Applying  the  above  theorem  to  equations 
(1.2)  and  (1.1),  we  conclude  that  samples  of  a  two-dimensional  BLP  signal  on  a  non-uniform 
rectangular  grid  shown  in  figure  (2.3)  can  be  used  to  reconstruct  it  uniquely. 

2.2  Results  in  Bivariate  Polynomial  Interpolation  Theory 

As  we  saw  in  Theorem  (2.1)  of  the  previous  section,  one  way  to  interpolate  in  Iln,  the 

space  of  polynomials  of  total  degree  n,  is  to  choose  n  +  1  distinct  lines  /o,  ...,/n,  and  select  the 
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Figure  2.3:  Geometric  distribution  of  the  sampling  points  of 
Theorem  (2.2)  for  n  =  2  and  m  —  1  . 

interpolation  points  from  these  lines.  More  specifically,  if  we  choose  t  +  1  points  on  /,,  the  ith 
line,  then  the  interpolation  problem  is  guaranteed  to  have  a  unique  solution.  Since  in  most  signal 
processing  applications,  the  support,  regions  of  the  Fourier  coefficients  of  the  two-dimensional 
signals  under  consideration  are  rectangular,  we  are  primarily  interested  in  interpolation  in 
the  space  of  polynomials  whose  coefficients  have  rectangular  rather  than  triangular  region  of 
support.  Thus,  it  seems  natural  to  extend  the  theorems  of  the  previous  section  from  IIn,  the 
space  of  polynomials  of  maximum  degree  n  in  w  and  z,  to  !!(„_„),  the  space  of  polynomials 
of  maximum  degree  n  in  u>  and  maximum  degree  n  in  z.  As  we  will  see  in  this  section, 
generalization  of  Theorem  (2.1)  provides  us  with  a  spectrum  of  powerful  sampling  schemes 
for  two-dimensional  BLP  signals.  More  specifically,  the  extension  of  Theorem  (2.1)  in  section 


(2.2.1)  implies  that  if  we  choose  N  distinct  lines  of  slope  1  in  an  image  with  a  N  x  N  region  of 
support  in  the  Fourier  domain,  then  2 i  4-  1  arbitrary  samples  on  the  ith  line  are  sufficient  for 
reconstruction  of  the  image.  In  section  (2.2.2),  we  generalize  this  result  to  the  case  where  all  the 
sampling  lines  have  fixed  positive  integer  slopes.  Finally,  in  section  (2.2.3)  we  will  use  a  modified 
version  of  Bezout’s  theorem  to  generalize  the  results  of  section  (2.2.2)  to  the  case  where  our 
sampling  lines  have  different  rational  slopes.  Unlike  the  constructive  proofs  of  sections  (2.2.1) 
and  (2.2.2),  the  algebraic  geometric  approach  of  section  (2.2.3)  will  not  provide  an  algorithmic 
procedure  for  the  actual  reconstructions. 

The  theoretical  results  of  this  section  are  applied  to  the  problem  of  reconstruction  from 
threshold  and  sine-wave  crossings  and  a  variety  of  other  reconstruction  problems  in  section 
(2.3).  The  actual  reconstruction  algorithms  and  their  numerical  properties  are  discussed  at 
length  in  Chapters  3  and  5. 

2.2.1  Sampling  Signals  on  Lines  of  Slope  1 

In  this  section,  we  will  extend  Theorem  (2.1)  of  section  (2.1)  to  the  case  where  the  inter¬ 
polation  is  done  in  ri(nn),  the  space  of  polynomials  of  maximum  degree  n  in  x  and  maximum 
degree  n  in  y.  This  is  because,  in  most  signal  processing  applications,  we  are  usually  interested 
in  images  with  square  support  in  the  Fourier  domain.  While  Theorem  (2.1)  can  be  applied  to 
these  signals,  since  many  of  the  coefficients  which  correspond  to  a  triangular  support  are  zero, 
the  required  number  of  interpolation  points  exceeds  the  number  of  unknown  Fourier  coefficients. 

Our  first  result  on  interpolation  in  II(n  n)  enables  us  to  choose  the  interpolation  points  on 
lines  which  pass  through  the  origin.  It  can  be  stated  in  the  following  manner: 


Theorem  2.3  Let  Iq,  be  distinct  lines  with  the  ith  line,  defined  by 


z  =  a,u) 


Q,-  ?  0 


and  consider  arbitrary  distinct  points  on  /,  given  by 


{(“'j*1.  I  3  =  °>  -)2<}  (2-3) 

where  the  ordering  of  the  lines  is  arbitrary.  If  none  of  the  interpolation  points  is  equal  to  (0,0), 
the  common  intersection  of  all  lines,  then  for  any  data  set 

{tf  |  °  <j  <  2«;0<  i  <  n}  (2.4) 

there  is  a  unique  bivariate  polynomial  of  the  form 

p(w,z)  =  (2-5) 

i=0  j=0 


such  that 


*{°)  =  0  <  3  <  2i;  0  <  i  <  n 


Proof:  Substituting  the  equation  of  the  kth  line,  4  into  p(ui,  z)  we  get 

n  n 

p(u),  akw)  =  a(t,j)tn’(afcu))J 

isOj'sd 

2« 

=  £ 


where 


i<*'  = 


^  a(t  —  m,  m)a™  0  <  t  <  n 

m=0 

n 

y*  a(i  —  m,  m)a£*  n  <  t  <  2r 


For  an  arbitrary  integer  a  >  —  1,  we  can  split  the  summation  in  equation  (2.7)  in  the  following  manner: 

r.  =  p(u<, akw)  - £fcjfc,u>’  -  £  (2-9) 


i  =  2n  —  n 


Setting  a  =  -1,  using  the  2n  +  1  points  of  ln,  we  can  uniquely  determine  6[n|  for  0  <  t  <  2 n.  In 
particular,  considering  equation  (2.8),  the  values  of  ip"*  and  enable  us  to  determine  o(0, 0),  a(n,  n), 
and  for  0  <  j  <  n  —  1.  This  is  because,  from  equation  (2.8),  we  have: 


jl 


1 


$ 

; 

s 


Similarly,  by  setting  s  =  0  in  equation  (2.9)  and  using  the  2n-  1  points  on  i„-i,  we  can  uniquely  specify 
11  for  1  <  t  <  2n  —  1.  This  can  be  done  because  the  determinant 


xn 

2 

2n  —  1 

*0 

*  *  *0 

0 

2n  —  1 

X\ 

Xl 

.  .  xA 

*2n-2 

2 

2n—  1 

*2n-2  * 

•  *  *2n-2 

is  non-sero  as  long  as  the  z.’s  are  different  from  each  other  and  from  aero.  Now  we  can  utilise  the 
values  of  b[n)  and  b[n  11  together  with  equation  (2.8)  to  find  a(0, 1)  and  a(l,0).  More  specifically,  from 
equation  (2.8)  we  have: 

b[k)  =  afca(0,l)  +  a(l,  0)  (2.13) 

Letting  k  —  n,  n  —  1  in  the  above  equation,  we  can  uniquely  specify  a(0, 1),  a(l,  0)  and  hence  fcj  for 
0  <  k  <  n  —  2. 

In  a  similar  manner,  the  values  of  ijn-i  and  &2n-V  can  be  used  to  find  a(n,  n  —  1)  and  a(n  —  l,n) 

(  t.  i 

and  hence  b2n_  x  More  specifically  from  equation  (2.8)  we  get: 

b2kn-i  =  a£a(n-l,n)  +  *a(n,  n  -  1)  (2.14) 

The  determinant  of  the  above  system  of  equations  for  k  =  n,  n  —  1  is  given  by 

a"-l  an-l  _  an-l  1  (o  l 

a"  a”-1  ”  1  (215) 

Taking  into  account  that  a*  #  0,  the  above  determinant  is  guaranteed  to  be  non  zero.  Thus  the 
coefficients  a(n,  n  -  1)  and  a(n  -  1,  n)  can  be  specified  uniquely. 

Repeating  the  above  procedure  for  s  =  l,...,2n  -  1,  we  can  find  all  the  coefficients  a(t ,j).  More 
specifically,  at  the  sth  stage,  we  know  b\ kJ_ll  and  hi*1  for  0  <  k  <  n  and  using  the  2(n  -  s)  —  1  points  of 

the  line  we  can  find  ij"-"  11  for  s  +  1  <  i  <  2n  —  s  —  1.  These  values  will  enable  us  to  uniquely 

specify  the  coefficients 

{o(t,j)|t  +  j  =  s  +  l,2n-s-l) 

and  hence  h^n-n-i  an<l  ^l+'i  f°r  ®Sk<n  —  s  —  2.  Consequently  we  can  completely  and  uniquely 
determine  all  the  coefficients  of  p(w,  z).  □ 

Thus,  our  proof  not  only  shows  that  the  appropriate  set  of  interpolation  points  on  lines 
passing  through  the  origin  results  in  a  unique  solution,  but  also  provides  us  with  a  recursive 
method  to  find  the  coefficients.  The  geometrical  distribution  of  the  interpolation  points  required 
by  this  theorem  for  n  =  2  is  shown  in  figure  (2.4). 

We  will  now  use  Theorem  (2.3)  to  define  a  sampling  strategy  for  BLP  signals  of  the  form 
given  by  equation  (1.1).  Since  g(WJt  W2)  of  equation  (1.2)  is  a  bivariate  polynomial  of  maximum 
degree  2 N  in  Wj  and  maximum  degree  2 N  in  W2,  invoking  Theorem  (2.3),  we  can  reconstruct 


27 


Figure  2.4:  Geometric  distribution  of  the  sampling  points  of 
Theorem  (2.3)  for  n  —  2. 

it  by  choosing  our  sampling  points  along  2 N  +  1  distinct  lines  given  by: 

W2  =  a,W!  0  <  i  <2N  (2.16) 

where 

Wx  = 

W2  =  e,2xy 

If  we  let 

a  = 


then  the  lines  in  the  -  W2  plane,  given  by  equation  (2.16),  will  correspond  to  lines  of  the 


in  the  x  -  y  plane.  Putting  all  of  this  together,  we  get: 


Corollary  2.1  Consider  a  bandlimited,  continuous  time,  periodic  signal  f(x,  y)  with  period  one 
in  the  x  and  y  directions  and  Fourier  series  representation 


N  N 


/(*.»)  = 

Y  E  ^2'v"2 

(2.18) 

Let  /o ,  -Jin  be  distinct  lines  in 

nl=-N  n2=-N 

the  x  -  y  plane  with  the  ith  line,  given  by 

y  =  x  +  ft 

(2.19) 

and  the  set  of  arbitrary  distinct 

samples  on  /,  given  by 

{(i',),y‘,))  1  o  <  >  <  2*'} 

(2.20) 

Then,  for  any  data  set 

{ 

t(p  j  0  <  j  <  2t  ;  0  <  t'  <  2JV} 

(2.21) 

we  can  reconstruct  f(x,y)  uniquely. 


An  example  of  the  geometric  distribution  of  the  sampling  points  required  by  this  theorem 
for  N  =  1  is  shown  in  figure  (2.5).  Figure  (2.5a)  shows  three  sampling  lines  of  unit  slope  in 
the  periodically  extended  version  of  a  BLP  signal,  and  figure  (2.5b)  shows  how  sampling  lines 
“wrap  around”  one  period  of  the  signal.  From  Corollary  (2.1)  we  can  conclude  that  the  set  of 
points  on  lines  of  slope  one  given  by : 

{(a^yf)  I  yf  =  xf  +  /?(,),  0  <  <  2,  o  <  j  <  2  «•} 

are  sufficient  to  guarantee  the  unique  reconstruction  of  the  signal. 

As  we  will  see  in  section  (2.3),  Corollary  (2.1)  can  be  applied  to  the  problem  of  reconstruction 
of  multidimensional  signals  from  threshold  crossings  by  choosing  the  interpolation  points  to  be 
the  intersections  of  level  crossings  and  the  sampling  lines. 


2.2.2  Sampling  Signals  on  Lines  of  Positive  Integral  Slope 

In  this  section,  we  will  generalize  Theorem  (2.3)  to  the  case  where  our  interpolation  points 
awe  chosen  on  curves  of  the  form  z  =  awm  as  opposed  to  straight  lines  passing  through  origin. 
Thus  the  results  in  the  previous  section  become  a  special  case  of  the  ones  we  will  derive  in  this 
section.  A  generalized  version  of  Theorem  (2.3),  can  be  stated  in  the  following  manner: 

Theorem  2.4  Let  cq,...,cp  be  distinct  curves  with  c,-,  the  ith  curve  given  by 

z  =  a,u;m  a,  ^  0 

where  m  <  n  is  an  arbitrary  positive  integer,  p  is  an  integer  satisfying 

p 

£2[(m  +  l)n  -  2mt  +  1]  >  (n  +  l)2 

«=o 

and 

=  0,...,(m+  l)n  -  2mt;  }  (2.22) 

is  a  set  of  distinct  points  on  li.  If  none  of  the  interpolation  points  defined  by  (2.22)  are  equal 
to  (0,0),  the  common  intersection  of  all  the  curves,  then  for  any  data  set 

2mi;  t'  =  0,  ...,p;} 

there  is  a  unique  bivariate  polynomial  of  the  form 

p(w,z)  =  ^^a(i',j')w‘zJ 
»=0;=0 

such  that 

p(ttfj*,,rj*))  =  tV  j  =  0,...,(m  +  l)n-2mt  t=0,...,n 

The  proof  is  included  in  Appendix  (A).  The  proof  not  only  shows  that  interpolation  points 
on  curves  of  the  form  z  =  a wm  result  in  a  unique  solution,  but  also  provides  a  recursive 
method  to  find  the  coefficients.  It  is  worthwhile  to  mention  that  unlike  Theorem  (2.3),  the 
number  of  interpolation  points  used  by  Theorem  (2.4)  might  be  larger  than  the  number  of 
unknown  coefficients.  More  specifically,  it  can  be  shown  that  the  number  of  sampling  curves, 
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Figure  2.6:  Geometric  distribution  of  the  sampling  points  of 
Theorem  (2.4)  for  n  —  2  and  m  —  2  . 


p  +  1,  is  given  by 


,n  +  1.. 

p  +  1  =  r— i 

m 


where  \k]  is  defined  to  be  the  smallest  integer  equal  to  or  larger  than  k.  Thus,  the  number 
of  interpolation  points  is  equal  to  the  number  of  coefficients  only  when  n  +  1  is  divisiole  by 
m.  The  geometrical  distribution  of  the  interpolation  points  required  by  the  above  theorem  for 
n  =  2,  m  =  2  is  shown  in  figure  (2.6). 

We  will  now  use  Theorem  (2.4)  to  define  a  generalized  sampling  strategy  for  BLP  signals  of 
the  form  given  by  equation  (1.1).  Since  g(W\,  W2)  in  equation  (1.2)  is  a  bivariate  polynomial  of 
maximum  degree  2 N  in  and  maximum  degree  2 N  in  W 2,  invoking  Theorem  (2.4),  we  can 
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reconstruct  it  by  choosing  our  sampling  points  along  2 N  +  1  distinct  lines  given  by 


W2  =  OiWf  0  <  »  <  2N  (2.23) 

where 

IV,  =  e*'2” 

W2  =  e>2*v 

If  we  let 

a  =  e>2'* 

then  the  lines  in  the  W\  -  plane  given  by  equation  (2.16)  will  correspond  to  lines  of  the 
form 

y  =  mx  +  fa  (2-24) 

in  the  i  -  y  plane.  Putting  all  of  this  together,  we  get: 


Corollary  2.2  Consider  a  bandlimited,  continuous  time,  periodic  signal  f(x,y)  with  period  one 
in  the  x  and  y  coordinates  and  Fourier  series  representation 

/(*,»)  =  E  E  F{n  1,n2)^2”n*  (2.25) 

nl=~N  n2= -N 

Let  Iq,  ...,lp  be  distinct  lines  in  the  x  —  y  plane  with  /„  the  ith  line,  given  by 

y  =  mi  +  ft  (2.26) 

where  m  <  2 N  is  an  arbitrary  positive  integer,  p  is  an  integer  satisfying 

p 

£[2JV(m+l)  ~  2m«  +  1]  >  (21V  +  l)2 

i=0 


and  the  set  of  arbitrary  distinct  samples  on  li  is  given  by 

{(x}0,»?))  I  0  <  j  <  2(m  +  1  )N  -  2 mi} 

Then,  for  any  data  set 

Oj')  |  0  <  j  <  2 (m  +  l)N  -  2mi  ;  0  <  *  <  p) 
we  can  reconstruct  /( x,y)  uniquely. 


(2.27) 


(2.28) 


Two  examples  of  the  geometric  distribution  of  the  sampling  points  required  by  this  result 
for  N  =  1  are  shown  in  figure  (2.7).  In  the  first  example,  the  slope  of  the  sampling  lines  is  1, 
and  the  distribution  of  the  points  is  identical  to  that  of  Corollary  (2.1).  In  the  second  example, 
the  slope  of  the  sampling  lines  is  2,  and  the  number  of  interpolation  points  exceeds  the  number 
of  Fourier  coefficients.  In  general,  the  distribution  of  samples  required  by  Corollary  (2.2)  varies 
as  a  function  of  the  slope  of  the  sampling  lines,  providing  us  with  a  spectrum  of  sampling 
techniques.  For  instance,  if  our  signal  has  5x5  region  of  support  in  the  Fourier  domain,  any 
one  of  the  following  sampling  sets  can  be  used  to  reconstruct  it  uniquely. 

1.  The  set  of  25  points  on  lines  of  slope  one  given  by 

{(x^1,  y|l))  |  =  xf  +  0  <  I  <  5,  0  <  j  <  4N  +  1  -  2i} 

2.  The  set  of  29  points  on  lines  of  slope  two  given  by 

{(*S‘).y]‘))  i  vf  =  2  xf  +  £(<),  0  <  «  <  2,  0  <  j  <  6AT  -  4.  +  1} 

3.  The  set  of  28  points  on  lines  of  slope  three  given  by 

{(x^.yj0)  |  yj°  =  3  xf  +  /?(*>,  0  <  »  <  1,  0  <  j  <  8JV  -  6t  +  1} 

4.  The  set  of  34  points  on  lines  of  slope  four  given  by 

{(x'^yj0)  |  yf  =  4  xf  +  /?<*>,  0  <  i  <  1,  0  <  j  <  ION  -  8,  +  1} 

5.  The  set  of  25  points  on  a  line  of  slope  5  given  by 


Note  that  for  cases  (1)  and  (5),  the  number  of  unknown  Fourier  coefficients  is  equal  to  the 
number  of  interpolation  points  whereas,  for  cases  (2),  (3)  and  (4),  we  need  more  samples  than 
coefficients.  As  we  will  see  in  section  (2.3),  the  results  derived  in  this  section  can  be  applied 
to  the  problem  of  reconstruction  from  level  crossings  by  choosing  the  intersection  of  sampling 
lines  and  threshold  contours  as  our  interpolation  points. 

2.2.3  Sampling  Signals  on  Lines  of  Rational  Slope 

In  the  previous  two  sections,  we  found  sufficient  conditions  for  reconstruction  of  multidi¬ 
mensional  BLP  signals  from  their  samples  on  lines  with  fixed  positive  integer  slopes.  In  this 
section,  we  use  a  modified  version  of  Bezout’s  theorem  to  derive  a  more  general  result  on  bi¬ 
variate  interpolation.  This  result  is  then  applied  to  find  the  sufficient  conditions  for  sampling 
BLP  signals  along  lines  with  different  rational  slopes.  Furthermore,  our  algebraic  geometric  ap¬ 
proach  enables  us  to  determine  the  sampling  conditions  under  which  the  reconstruction  problem 
is  guaranteed  to  have  infinitely  many  solutions. 

As  we  will  see,  the  theoretical  results  developed  in  this  section  are  more  general  than  the 
ones  in  the  previous  sections.  In  fact,  one  can  show  that  all  the  theorems  of  sections  (2.2.1) 
and  (2.2.2)  are  special  cases  of  the  results  in  this  section.  However,  unlike  the  constructive 
proofs  of  the  previous  sections,  the  algebraic  geome^-ic  approach  does  not  provide  us  with  an 
algorithmic  procedure  for  the  actual  reconstructions. 

Let  us  begin  with  the  important  concept  of  irreducibility.  A  bivariate  polynomial  is  said  to 
be  irreducible  over  complex  numbers,  if  it  cannot  be  factored  into  polynomials  of  smaller  degree 
with  complex  coefficients.  Thus,  two  irreducible  polynomials  have  no  common  factors  unless  the 


one  with  smaller  total  degree  is  a  factor  of  the  other  one.  Bezout’s  theorem  is  concerned  with 


the  number  of  common  zeros  of  two  bivariate  polynomials  and  can  be  stated  in  the  following 
manner: 

Theorem  2.5  (Bezout  [37,38])  If  two  bivariate  polynomials  of  total  degree  r  and  s  given  by: 

p(x,y)  = 

t=0  j=o 

<iix,y)  = 

«=0 j=0 

have  no  common  factors  of  degree  greater  than  zero,  then  they  have  at  most  rs  common  zeros. 

In  this  theorem,  the  total  degree  of  a  polynomial  in  two  variables  is  defined  in  terms  of 
the  sum  of  the  degrees  in  each  variable.  That  is,  the  total  degree  of  the  two-dimensional 
polynomial  p(x,  y)  is  equivalent  to  the  degree  of  the  one-dimensional  polynomial  p(x,x).  Since 
we  are  primarily  interested  in  signals  with  rectangular  Fourier  domain  support,  we  need  to 
modify  Bezout’s  theorem  in  order  to  find  a  tighter  upper  bound  on  the  number  of  common 
zeros  of  two  polynomials  whose  coefficients  have  rectangular  regions  of  support.  The  modified 
Bezout’s  theorem  originally  derived  by  Zakhor  and  Izraelivitz  [39]  can  be  stated  in  the  following 
way: 

Theorem  2.6  (Modified  Bezout’s  Theorem  [39])  Consider  two  bivariate  polynomials  p(x,  y)  € 
n(Ar,,Ar„)  and  q(x,y)  €  n(Mi My)  of  the  form 

N,  ”, 

pix>y)  =  £5Za(*J)*V 
1=0 >=0 

M,  M, 

i(x>y)  =  I 

«=0j=0 

If  p  and  q  have  no  common  factors  of  degree  greater  than  zero,  then  they  have  at  most  NzMy  +  MXNV 
common  zeros. 


The  proof 
n(Af,  ,/v,)  with 
zeros  with  an 


winr  w- 
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is  included  in  Appendix  (B).  Theorem  (2.6)  implies  that  if  a  polynomial  p(x,  y)  € 
maximum  degrees  Nx  in  x  and  Nv  in  y,  has  more  than  NxMy  +  MXNV  common 
irreducible  polynomial  q(x,y)  €  with  maximum  degrees  Mx  in  x  and 


Mv  in  y,  then  q(x,  y)  must  be  a  factor  of  p( x,y).  This  consequence  of  Theorem  (2.6)  can  be 
used  to  derive  the  most  general  theorem  of  this  chapter: 


Theorem  2.7  Let  co,Ci,...,cp  be  distinct  bivariate  irreducible  polynomials  with  the  maximum 
degrees  of  c,  in  w  and  z  given  by  and  and  p  being  an  integer  satisfying  either  of  the 
following  two  conditions: 

p 

(2.29) 

t=0 

(2-3°) 

i=0 

Define  A,  to  be  the  set  of 

5(0  =  m^(nw  -  £  mW)  +  -  £  m<fc>)  +  1 

k=0  k=0 

points  on  c,  given  by 

A,  =  {(u^,  zf)  |  =  0,  0  <  j  <  S(0>  (2.31) 

If  none  of  the  interpolation  points  given  by  (2. SI)  are  on  the  intersections  of  two  or  more  of 
the  Ci ’s,  then  for  any  data  set 

{<;,)  i  o  <  *  <  P,  o  <  j  <  5(0} 

there  is  a  unique  bivariate  polynomial  of  the  form 

p(w,z)  = 

«=o  >=o 

such  that 

p(^\  zf)  =  tf  0  <  «  <  p,  0  <  j  <  5(0}  (2.32) 


Proof:  To  show  that  there  is  a  unique  polynomial  which  satisfies  (2.32),  we  have  to  show  that 

there  are  no  polynomials  in  j  which  vanish  at  all  the  interpolation  points  uAt.  Suppose,  on  the 

contrary,  that  there  is  a  polynomial  q(w,z)  €  which  vanishes  at  all  the  interpolation  points. 
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Since  g  has  mi?1***  4-  nwm[0)  +  1  common  zeros  with  Co,  by  the  modified  version  of  Bezout’s  theorem, 
Co  must  be  a  factor  of  q(w,z).  That  is 

<7(u>,2)  =  c0(u>,z)gm(ti;,  z) 

where  g11*^,  z)  is  a  polynomial  of  maximum  degree  n„,  -  m??1  in  w  and  nz  -  mi"1  in  z.  Furthermore, 
since  by  hypothesis,  none  of  the  interpolation  points  on  Cj  are  on  c0  and  q(u/,z)  has  1  +  — 

m«? 1 )  +  mi, 1  (nz  —  m*  )  common  zeros  with  ci,  it  must  be  true  that  the  same  number  of  common 
zeros  with  c j.  Taking  into  account  the  irreducibility  of  cj,  by  modified  version  of  Bezout’s  theorem,  c* 
must  be  a  factor  of  g(1|(to,  z )  and  hence  g(«n,  z)  . 

Repeating  the  above  argument  for  C2,  we  get: 

q(ui,z)  =  c„(tu,z)  ...  cv.i(w,z)  q[,,i(w,z) 

i'- 1  i»-i 

where  g(,,,(tn,  z)  has  maximum  degree  nw  ~  in  w  and  maximum  degree  n2  —  in  z  and 

«=<>  >=o 

»•- 1  /-- 1 

has  1  +  m?k|)  +  mi,’,(nv,  -  ^  mif1)  common  zeros  with  c,,.  Since  c,,  is  irreducible,  by 

*=<l  k=l) 

modified  version  of  Bezout’s  theorem,  it  must  be  a  factor  of  glpl(u>,  z).  This  contradicts  our  hypothesis 
since  by  inequalities  (2.32)  and  (2.29),  the  degree  of  c|(  in  either  w  or  z,  is  larger  than  that  of  gtp)(ti>,  z). 

□ 

The  above  theorem  is  the  generalized  version  of  Theorems  (2.3)  and  (2.4).  More  specifically, 
if  c,  is  an  irreducible  polynomial  of  the  form 

z  —  <X{Wl  (2.33) 

then  Theorem  (2.7)  simply  reduces  to  Theorem  (2.4)  for  positive  integer  values  of  l,  and  to 
Theorem  (2.3)  for  1=1.  This  is  because  polynomials  of  the  form  given  by  equation  (2.33) 
are  known  to  be  irreducible.  An  example  of  the  distribution  of  sampling  points  required  by 
Theorem  (2.7)  for  nw  =  nz  =  2,  mffl  =  =  1  and  =  2  is  shown  in  figure  (2.8). 

There  are  two  classes  of  irreducible  polynomials  which  are  particularly  useful  in  deriving 
sampling  strategies  for  multidimensional  periodic  signals  of  the  form  given  by  equation  (1.1). 
These  two  classes  of  polynomials  in  W i  and  W?  are  of  the  form: 

W 2M»  =  a W™’  ,  Mz>  0,  Mv>  0  (2.34) 
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w“*w™*  =  a  ,  Mx  >  0,  M„  >  0 


(2*5) 


where  M„  and  Mx  are  positive  integers  which  are  relatively  prime  with  respect  to  each  other. 


Using  the  fact  that  W\  and  W2  of  equation  (1.1)  are  related  to  i  and  y,  the  signal  coordinates, 


-  J 


Wi  = 


W2  =  e>2*v 


and  letting 


a  = 


the  curves  in  the  W j  -  W2  plane  given  by  equations  (2.34)  and  (2.35)  correspond  to  lines  of  the 


A/„y  =  0  +  Mzx  ,  Mx,  My  >  0 


(2.36) 


Mvy  +  Mz x  =  0  ,  Mx,  My  >  0 


(2.37) 


in  the  x  -  y  plane.  Geometrically,  equations  (2.36)  and  (2.37)  correspond  to  lines  with  positive 


and  negative  rational  slopes  respectively.  Therefore,  lines  with  positive  and  negative  rational 


slopes  in  the  x-y  domain  correspond  to  irreducible  bivariate  polynomials  in  the  Wj  -  W2  plane. 


We  can  use  this  fact,  together  with  Theorem  (2.7)  to  define  the  following  sampling  strategy: 


Corollary  2.3  Consider  a  bandlimited,  continuous  time,  periodic  signal  f{x,y)  with  period  one 
in  x  and  y  direction  and  Fourier  series  representation: 


N  N 

/(*,»)  =  £  E  nil,**)**2™1  ej2rvn2 

nl  =  -N  n2=-N 


(2.38) 


Let  Iq,  ...,/p  be  distinct  lines  in  the  x-y  plane  with  li,  the  ith  curve  given  by 


M{']y  =  M&x  +  pi 


(2.39) 
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Figure  2.8:  Geometrical  distribution  of  the  sampling  points 
of  Theorem  (2.7)  for  nw  -  nz  —  2. 


where  .  nd  are  positive  or  negative  integers  which  are  relatively  prime  with  respect  to 
each  other.  Let  p  be  an  integer  satisfying  either  one  of  the  following: 


2*<£K,| 

»=o 

(2.40) 

2N<'£\mM\ 

i=0 

(2.41) 

Suppose  that  the  set  of 

S(i)  =  Ja/W 

|(21V-2K*)|)  +  |AfW|(2JV-2Kk)|)  +  1 

*=o  *=o 

arbitrary  distinct  samples  on  li,  is  given  by 

1  o  <  »'  <  P,  o  <  j  <  5(»)} 

(2.42) 

If  none  of  the  interpolation  points  given  by  equation  (2. 42)  are  on  the  intersection  of  two  or 
more  of  the  c,  ’s,  then  for  any  data  set 


{tf  I  0  <  j  <  5(.)}  (2.43) 
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Figure  2.9:  Geometric  distribution  of  the  sampling  points  of 
Corollary  (2.3)  for  N  —  1. 

we  can  reconstruct  f(x,y)  uniquely. 

An  example  of  the  geometric  distribution  of  the  sampling  points  required  by  this  theorem 
for  N  =  1  and  M^0'  =  -1,  My0'  =  1  ,  M^'  =  1,  My1'  =  2  is  shown  in  figure  (2.9).  The  above 
result  is  more  general  than  our  previous  results  in  the  sense  that  Corollary  (2.1)  becomes  a 
special  case  of  it  when 

mJ*'  =  =  1  ,  vi 

and  Corollary  (2.2)  becomes  its  special  case  when 


All  the  theorems  we  have  discussed  so  far,  provide  us  with  sufficient  conditions  for  unique 
recovery  of  polynomials  or  their  associated  BLP  signals.  We  can  use  arguments  similar  to  the 
proof  of  Theorem  (2.7)  to  find  conditions  under  which  the  reconstruction  problem  will  have 
infinitely  many  solutions: 

Theorem  2.8  Let  co,ci,...,cp  be  distinct  bivariate  polynomials  with  the  maximum  degrees  of  Cj 
in  w  and  z  given  by  and  where  p  is  an  integer  satisfying  the  following  two  inequalities: 

nw  >  m»)  (2-44) 

t=0 

n,  >  53  mi*'  (2.45) 

t=0 

Suppose  that  we  choose  S(i),  an  arbitrary  number  of  interpolation  points,  on  the  ith  curve,  Cj. 
Then  for  any  data  set 

I  0  <  «  <  p,  0  <  j  <  S(i)} 

there  are  an  infinite  number  of  bivariate  polynomials  of  the  form 

p(w,z)  =  53  53a(.\j)ti/V 
»=o  >= o 

such  that 

=  ty  ’  0  <  i  <  p,  0  <  j  <  S{i)  (2.46) 

Proof:  The  proof  is  similar  to  the  proof  of  Theorem  (2.7).  To  show  that  there  are  an  infinite 

number  of  polynomials  which  satisfy  (2.46),  it  is  sufficient  to  show  that  there  is  at  least  one  polynomial 
of  maximum  degree  <  n„,  in  w  and  <  nz  in  z  which  vanishes  at  all  the  interpolation  points  defined  by 
the  theorem.  The  most  obvious  choice  for  this  polynomial  is 

p 

g(w>z)  =  n  *<•.«) 

1=0 

Since  by  inequalities  (2.44)  and  (2.45),  the  maximum  degree  of  q(w,  z)  is  less  than  nw  in  w  and  less  than 
n,  in  2,  there  exists  infinite  number  of  polynomials  in  n(„w.nj|  which  satisfy  (2.46).  □ 

TVanalating  the  above  result  to  the  signal  domain  via  our  usual  transformations: 


and  assuming  the  irreducible  curves  of  the  above  corollary  to  be  of  the  form  shown  in  equations 
(2.34)  and  (2.35),  we  arrive  at  the  following  result: 


Corollary  2.4  Consider  a  bandlimited,  continuous  time,  periodic  signal  f{x,y)  with  period  one 
in  x  and  y  direction  and  Fourier  series  representation: 

f(x,y)  =  Z  £  F(ni,n2)^^  ej2lry"2  (2.47) 

nl=-iV  n2=-N 

Let  lo,  be  distinct  lines  in  the  x  -  y  plane  with  L,  the  ith  curve  given  by 

M^y  =  M|‘>i  +  0i  (2.48) 

where  Mi''1  and  are  positive  or  negative  integers  which  are  relatively  prime  with  respect  to 
each  other.  Let  p  be  an  integer  satisfying  the  below  two  inequalities: 

2JV>£|m<<)|  (2.49) 

i=0 

2N>'£\mM\  (2.50) 

«=o 

Suppose  that  we  choose  S(i),  an  arbitrary  number  of  interpolation  points  on  the  ith  line  /j. 
Then,  for  any  data  set 

(4°  I  0  <  j  <  S{i)}  (2.51) 

there  are  infinitely  many  functions  of  the  form  given  by  equation  (2. 4  7)  such  that: 

/(* ■J’),  yj‘))  =  4°  0  <  i  <  p,  0  <j<  S(i) 

In  terms  of  sampling  multidimensional  periodic  signal,  the  above  theorem  implies  that 
unless  the  number  of  interpolation  lines  and  their  respective  slopes  are  chosen  carefully,  we 
might  encounter  situations  where  the  interpolation  problem  has  non  unique  solutions. 

To  summarize  this  section,  we  have  derived  a  variety  of  results  on  multivariate  polynomial 
interpolation  in  n(n  n).  Exploiting  the  polynomial  representation  of  multidimensional  BLP 
signals,  we  applied  these  results  to  the  problem  of  reconstruction  from  non-uniformly  spaced 


samples.  More  specifically,  our  results  provide  us  with  sufficient  conditions  for  the  unique 
reconstruction  of  signals  from  their  samples  on  various  lines  of  positive  or  negative  rational 
slopes.  A  summary  of  all  the  theorems  and  corollaries  of  the  last  two  sections  is  included 
at  the  end  of  this  chapter.  In  the  next  section,  we  will  apply  these  results  to  a  variety  of 
multidimensional  reconstruction  problems. 

2.3  Applications  of  the  Line-Sampling  Strategy 

In  section  2.2,  we  found  sufficient  conditions  under  which  samples  of  a  multidimensional 
signal  can  be  used  to  uniquely  specify  it.  More  specifically,  the  theoretical  results  of  the  previ¬ 
ous  section  provide  us  with  a  variety  of  distributions  of  sampling  points  along  lines  of  positive 
or  negative  rational  slope  which  result  in  unique  recovery  of  a  signal.  Although  our  main 
motivation  for  deriving  these  semi-implicit  sampling  strategies  has  been  the  problem  of  recon¬ 
struction  from  multiple  level  threshold  crossings,  they  can  also  be  applied  to  a  variety  of  other 
reconstruction  problems  with  non-uniformly  spaced  samples  in  areas  such  as  machine  vision 
[18,19],  radio  astronomy  [20]  and  computed  tomography  [21],  as  well  as  natural  sciences  such 
as  geology,  meteorology,  and  oceanography.  These  problems  can  be  divided  into  the  following 
three  categories: 

1.  Recovery  of  signals  from  their  crossings  with  arbitrary  functions.  A  special  case  of  this 
problem  is  reconstruction  of  signals  from  their  single  or  multiple  level  threshold  crossings. 
Another  special  case  is  reconstruction  of  signals  from  their  crossings  with  periodic  func¬ 
tions,  which  has  potential  applications  in  the  conversion  of  halftone  images  to  continuous 
tone  ones. 

45 


2.  Recovery  of  signals  from  their  projections.  Applications  of  this  problem  include  such 


diverse  fields  as  X-ray  tomography,  transmission  electron  microscopy  and  radio  astronomy. 

3.  Recovery  of  signals  whose  values  are  known  along  certain  paths,  curves  or  contours. 
This  problem  arises  frequently  in  the  area  of  machine  vision,  where  depth  or  distance 
information  is  available  on  the  zero  crossing  contour  of  the  convolution  of  the  image  with 
the  Laplacian  of  a  Gaussian. 

Although  the  main  focus  of  this  thesis  is  the  problem  of  reconstruction  from  multiple  level 
threshold  crossings,  in  this  section,  we  will  briefly  describe  the  way  our  semi-implicit  schemes 
of  section  (2.2)  can  be  applied  to  some  of  these  other  problems. 

2.3.1  Recovery  from  Crossings  with  Arbitrary  Functions 

In  reconstructing  signals  from  their  crossings  with  arbitrary  functions,  the  intersections  of 
sampling  lines  with  the  crossing  contours  are  used  as  interpolation  points.  The  steps  involved 
in  the  reconstruction  of  BLP  signals  from  their  multiple  level  crossings  can  be  summarized  as 
follows: 

1.  Find  the  level  crossing  contours  associated  with  the  thresholds. 

2.  Find  a  set  of  sampling  lines  with  rational  slopes  whose  intersections  with  the  crossing 
contours  satisfies  the  distributions  required  by  Corollaries  (2.1),  (2.2)  or  (2.3). 

3.  Use  the  intersections  to  find  the  coefficients  of  the  polynomial  associated  with  the  signal 
or  equivalently  the  Fourier  coefficients  of  the  signal. 

In  Chapter  3,  we  will  propose  a  variety  of  techniques  for  carrying  out  the  third  step  of  the  above 


difficulty  in  applying  the  above  procedure  to  actual  reconstruction  problems  is  that,  in  general, 
we  are  not  guaranteed  to  get  enough  intersections  between  sampling  lines  and  threshold  contours 
to  satisfy  the  distribution  requirements  of  Corollaries  (2.1),  (2.2)  or  (2.3).  Of  course,  as  will 
be  seen  in  section  (3.4),  we  can  find  guidelines  which  help  us  choose  the  slope  and  position 
of  our  sampling  lines  in  an  “optimal”  fashion  so  that  the  number  of  intersections  of  sampling 
lines  with  the  threshold  contours  is  maximized.  However,  the  problem  still  remains  that  for 
small  number  of  thresholds,  nt,  we  might  not  be  able  to  find  any  set  of  sampling  lines  which 
satisfy  our  theoretical  requirements.  Indeed,  this  becomes  our  main  motivation  for  deriving  the 
less  restrictive  results  of  Chapter  4.  Meanwhile,  we  can  make  a  few  observations  for  the  more 
general  problem  of  reconstruction  from  crossings  with  arbitrary  functions.  As  it  turns  out,  for 
certain  class  of  functions  such  as  sinusoids,  the  number  of  intersections  of  sampling  lines  with 
function  crossings  of  the  signal  becomes  signal  independent.  As  an  example,  consider  the  eye 
picture  shown  in  figure  (2.10)  and  its  sinusoid  crossings  with  the  function: 

h(x,y)  =  A(l+  cos{2n(px  +  qy)))  (2.52) 

shown  in  figure  (2.11).  As  seen,  the  above  sinusoid  assumes  its  maximum  value  2 A  and  minimum 
value  0  on  equidistant  lines  of  slope  - 1 .  Thus,  if  we  impose  the  condition 

0  =  [h(x,  y)]min  <  f{x,y)  <  IMZ>  y)]mai  =  2A  ,  V(x,  y)  (2.53) 

on  the  amplitude  of  our  signal,  then  between  every  two  neighboring  parallel  lines  where  h(x,  y) 
assumes  its  minimum  and  maximum  values,  there  exists  a  contour  of  crossings  of  h(x,  y)  and 
our  signal  f(x,y).  Since,  the  intensity  of  the  eye  picture  lies  in  the  range  [0,256],  to  satisfy 
inequality  (2.53),  the  value  of  A  for  the  crossings  shown  in  figure  (2.11)  was  chosen  to  be  128. 
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If  we  now  sample  these  sinusoid  crossings  of  our  signal  along  lines  of  rational  slope  ^  given  by 


l  :  y  =  —i  + 
m 


we  are  guaranteed  to  get  exactly  2[pm  +  qn)  samples.  As  shown  in  figure  (2.11),  each  sample 
on  line  /  is  located  between  the  two  intersections  of  the  sampling  line  with  two  lines  at  which 
h(x,  y)  assumes  its  maximum  and  minimum  values. 

We  can  use  the  above  example,  together  with  the  distribution  requirements  of  Corollaries 
(2.1),  (2.2)  and  (2.3),  to  find  the  exact  number  of  sinusoids  needed  for  recovery  of  a  signal  via  a 
given  set  of  sampling  lines.  For  instance,  suppose  that  we  are  interested  in  sampling  crossings 
of  a  sinusoid  with  a(2JV  +  l)x(2iV  +  l)  BLP  signal  along  lines  of  rational  slope  ^  where  n  and 
m  are  relatively  prime  with  respect  to  each  other.  Invoking  Theorem  (2.7)  we  realize  that  for 
this  sampling  strategy  the  maximum  number  of  samples  needed  on  any  of  the  sampling  lines  is 
2 N(n  +  m)  +  1.  Since  the  number  of  intersections  of  a  sinusoid  with  period  i  in  the  x  direction 
and  ^  in  the  y  direction  with  a  line  of  slope  *  is  2  [pm  +  yn],  in  order  to  satisfy  the  sampling 
requirements  of  Theorem  (2.7),  we  must  have 

2(pm  +  qn)  >  2 N(m  +  n)  +  1 


(pm  +  qn)  >  (m  +  n)N  (2.55) 

A  necessary  condition  for  satisfying  the  above  inequality  is  that  at  least  one  of  the  quantities 
p  or  q  must  be  greater  than  N.  This  in  turn  implies  that  using  the  semi-implicit  sampling 
strategies  of  this  chapter  a  BLP  signal  can  be  uniquely  reconstructed  from  its  crossings  with 
a  sinusoid  only  if  the  frequency  of  the  sinusoid  lies  outside  the  bandwidth  of  the  signal.  More 
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Figure  2.10:  The  original  eye  picture  with  Zl  x  31  region  of  support  in  the  Fourier  domain. 


Maximum  lines 
Minimum  lines 


Figure  2.11:  Sinusoid  crossings  of  the  original  eye  picture.  The  lines  at  which  the  sinusoid 
assumes  its  maximum  and  minimum  values  are  also  shown. 

49 


v> 


4 


£ 


« 

y. 


sf. 

I 


ft 


« 

$ 


u~a  u-avm  wire  u-v v-v v/v  1  \  rar\  V1  V"  P*.'  VI KT  *\.^  K1  ) 


'. -mm ran*  y  wk'’»vr>i' 


generally,  reconstruction  of  an  (2 N  +  1)  x  (2 N  +  1)  signal  from  samples  of  its  crossings  with  L 
sinusoids  along  lines  of  slope  j*-  is  possible  if  and  only  if: 


1.  The  maximum  and  minimum  values  of  each  sinusoid  is  larger  and  smaller  than  the  max¬ 
imum  and  minimum  values  of  the  signal  under  consideration. 


2.  The  periods  of  the  sinusoids  along  x  and  y  directions  given  by  —  and  —  satisfy: 


^2{Pim  +  Vi*)  >  ( m+n)N 


(2.56) 


An  important  practical  application  of  reconstruction  from  crossings  with  periodic  functions 
is  recovery  of  contone  images  from  halftone  images  [17].  The  halftone  process  has  been  used  for 
more  than  a  century  for  converting  continuous  tone  pictures  into  regular  patterns  of  black  and 
white  dots  which  can  then  be  printed.  The  size  of  each  dot  is  related  to  the  tone  at  that  same 
place  in  the  original  image  being  reproduced.  Mathematically  speaking,  the  halftone  version 
of  a  continuous  tone  image  can  be  obtained  by  comparing  the  value  of  the  signal  with  a  two- 
dimensional  periodic  function  and  producing  a  white  or  black  pixel  on  a  high  contrast  medium 
depending  on  whether  the  signal  value  is  higher  or  lower  than  that  of  the  periodic  function. 
Thus,  if  the  period  of  the  thresholding  function  is  a  sinusoid  satisfying  the  inequality  (2.55),  we 
can  be  guaranteed  unique  reconstruction  of  the  continuous  tone  image  by  sampling  the  halftone 
image  along  lines  of  slope  rr.  A  few  examples  of  such  reconstructions  are  shown  in  Chapter  3. 


2.3.2  Recovery  from  Projections 


Another  major  application  of  the  theoretical  results  of  section  (2.2)  is  in  the  area  of  recon¬ 
struction  of  multidimensional  signals  from  their  projections.  In  this  section,  we  will  show  that 


m 


m 

Iw 
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the  one- project  ion  theorem  due  to  Mersereau  and  Oppenheim  [45]  is  a  special  case  of  Theorem 
(2.4)  and  can  be  generalized  to  situations  with  an  arbitrary  number  of  projections. 


Consider  a  bandlimited  two-dimensional  signal  of  order  M  and  bandwidth  W  whose  Fourier 


transform  is  given  by  [45]: 


FWl ,-l)  =  tn  E  E  /<  W ,  W  (-1  .< <*)  (2-57) 


w !  .  ^  ^ 

m=0  n=0 


where 


bw{ui,uj2)  - 


1  lu^l  <  1,  |w2|  <  1 


(2.58) 


0  elsewhere 


The  key  result  in  recovering  signals  from  their  projections  is  the  projection  slice  theorem,  which 


essentially  relates  the  one-dimensional  Fourier  transform  of  a  projection  of  a  signal  to  a  slice  of 


its  two-dimensional  Fourier  transform.  In  many  applications,  the  only  available  information  for 


reconstructing  a  signal  is  samples  of  the  slices  of  its  two-dimensional  Fourier  transform.  The 


one-projection  theorem  of  Mersereau  and  Oppenheim  ([45])  basically  states  that  M2  samples  on 


one  slice  of  the  two-dimensional  transform  is  sufficient  for  unique  reconstruction  of  the  signal, 


provided  the  angle  of  the  projection  is  chosen  to  be  a  critical  angle.  To  prove  the  one-projection 


theorem  for  the  special  case  where  the  projection  angle  is  given  by  9  =  tan  1  M,  let 


Zl  =  e_; 


Z2  —  e 


-  „-lw“  2 


in  equation  (2.57).  Then,  the  Fourier  transform  shown  in  equation  (2.57)  can  be  written  as  a 


polynomial  in  terms  of  Z\  and  Z2 : 


2  M-1M-1 


nz,.z,)  = 


m=0  n=0 


R 


I 

& 


Thus,  invoking  Theorem  (2.4),  M2  samples  on  a  line  of  the  form: 


Z2  =  Z? 


in  the  Z\  -  Zi  domain,  or  equivalently 


<jJ2  =  MuJ\ 


in  the  —  u/2  plane  are  sufficient  for  unique  recovery  of  the  coefficients  of  F(o;],u;2).  Exploiting 
the  fact  that  Theorem  (2.4)  is  a  special  case  of  Theorem  (2.7),  we  can  generalize  this  result  to 
reconstruction  from  an  arbitrary  number  of  slices  with  rational  slopes.  More  specifically,  we 


Corollary  2.5  Consider  a  signal  whose  Fourier  transform  is  given  by  equation  (2.57).  Let 
So,  ...sp  be  distinct  slices  in  the  wj ,  w2  plane  with  Si  the  ith  slice  given  by: 

P2(,)w2  =  pj% 

where  P ^  an<^  Pi'^  are  positive  or  negative  integers  which  are  relatively  prime  with  respect  to 
each  other.  Let  L  be  the  smallest  integer  such  that  either  one  of  the  following  is  satisfied: 


M  -  1  <  £  \p2 


Suppose  that  the  set  of 


<3(0  =  !P2(,)l(M-i-LlA(fc)i)  +  I^KM-i-^lP^D  +  i 


arbitrary  distinct  samples  on  s,,  is  given  by 


{(^CO.^C/)  I  0  <  3  <  <3(0} 


Then,  if  none  of  the  interpolation  points  is  equal  to  (wi,w2)  =  (0,0),  coefficients  of  F(w i,w2), 
or  equivalently  our  original  signal  can  be  uniquely  reconstructed  from  the  samples  of  the  Fourier 
transform  of  the  signal. 


We  will  not  investigate  the  consequences  of  the  above  result  experimentally,  simply  because 
it  is  not  directly  related  to  the  main  focus  of  this  thesis,  which  is  reconstruction  from  multiple 
level  threshold  crossings. 

2.3.3  Recovery  from  Signal  Values  along  specific  Paths 

In  many  reconstruction  applications,  the  value  of  a  signal  is  known  along  certain  contours, 
paths  or  curves  and  it  is  necessary  to  recover  the  signal  from  this  information.  Clearly,  our 
results  of  previous  sections  can  be  applied  to  these  situations  by  using  the  intersections  of 
sampling  lines  and  the  contours  as  interpolation  points 

As  an  example,  consider  the  surface  interpolation  problem  which  arises  frequently  in  machine 
vision.  Computational  theories  of  structure  from  motion  [40]  and  stereo  vision  [41]  only  specify 
the  computation  of  three  dimensional  surface  information  at  special  points  in  the  image.  Yet, 
the  visual  perception  is  clearly  of  complete  surfaces.  To  account  for  this,  a  computational 
theory  for  interpolating  surfaces  from  visual  information  has  been  proposed  by  many  researchers 
[18].  More  specifically,  if  we  view  the  human  early  visual  system  as  a  symbolic  manipulator, 
we  can  consider  visual  processing  as  a  series  of  transformations  from  one  representation  to 
another  [42,43].  The  first  transformation  is  from  images  to  a  description  called  the  primal 
sketch  of  those  locations  at  which  the  image  irradiances  change.  These  locations  can  be  found 
by  finding  the  zero  crossings  of  the  convolution  of  the  original  image  with  the  Laplacian  of 
a  Gaussian.  Next,  primal  sketch  descriptions  of  several  images  are  matched,  either  by  the 
stereo  or  motion  computation,  to  obtain  a  description  of  the  surface  information  at  the  zero 
crossings.  This  representation  is  called  the  raw  2 sketch.  Finally,  the  raw  2 \D  sketch  is 
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interpolated  to  obtain  complete  surface  descriptions,  called  the  full  2  \D  sketch  [44].  Since  the 
input  representations  for  the  interpolation  step  consist  of  explicit  surface  information  such  as 
distance  or  relative  distance  along  the  zero  crossings  of  the  convolved  image,  the  intersections 
of  our  sampling  lines  with  the  zero  crossing  contours  of  the  convolved  image  can  be  ui^d  as 
interpolation  points  for  finding  the  polynomial  associated  with  the  surface  under  consideration. 
Under  these  circumstances,  the  theoretical  results  of  section  (2.2),  provide  us  with  a  variety  of 
sufficient  conditions  for  interpolating  the  distance  values  on  selected  points  of  the  zero  crossing 
contour  of  the  convolved  image. 

To  summarize  this  chapter,  we  have  developed  a  variety  of  results  on  interpolation  of  poly¬ 
nomials  in  n(n,n).  These  results  ultimately  have  been  used  to  derive  a  variety  of  semi-implicit 
sampling  strategies  for  reconstruction  of  multidimensional  signals  from  their  crossings  with  ar¬ 
bitrary  functions,  their  non-uniformly  spaced  samples,  and  their  projections.  A  summary  of 
all  the  theorems  and  corollaries  derived  in  this  chapter  are  included  in  table  (2.1).  In  the  next 
chapter,  we  will  propose  reconstruction  algorithms  for  the  specific  problem  of  reconstruction 
from  multiple  level  threshold  crossings. 


Theorem  or  Corollary 
Number 

Reference  for  proof 

Region  of 
Support 

Description 

Theorem  (2.1) 

Gasca  and  Maetzu  [36] 

triangular 

polynomial  interpolation 
with  n  +  1  lines  and  t  +  1 
samples  on  the  tth  line 

Theorem  (2.2) 

[29,32,33,34] 

rectangular 

polynomial  interpolation 
with  non-uniform 
rectangular  sampling 

Theorem  (2.3) 

section  (2.2.1) 

rectangular 

polynomial  interpolation 
with  n  +  1  lines  through  the 
origin  and  2t  +  1  samples 
on  the  tth  line 

Corollary  (2.1) 

Theorem  (2.3) 

rectangular 

recovery  of  BLP  signals 
from  their  samples  on 
lines  of  unit  slope 

Theorem  (2.4) 

appendix  (A) 

rectangular 

polynomial  interpolation 
from  samples  on  curves 
of  the  form  z  =  aw"‘ 

Corollary  (2.2) 

Theorem  (2.4) 

rectangular 

recovery  of  BLP  signals 
from  samples  on  lines 
of  positive  integer  slope 

Theorem  (2.5) 

Bezout  [37,38] 

triangular 

upper  bound  on  the  number 
of  common  zeros  of  two 
polynomials 

Theorem  (2.6) 

appendix  (B) 

rectangular 

upper  bound  on  the  number 
of  common  zeros  of  two 
polynomials 

Theorem  (2.7) 

section  (2.2.3) 

rectangular 

polynomial  interpolation  with 
samples  on  irreducible  curves 

Corollary  (2.3) 

Theorem  (2.7) 

rectangular 

recovery  of  BLP  signals 
from  their  samples  on 
lines  of  rational  slope 

Theorem  (2.8) 

section  (2.2.3) 

rectangular 

necessary  conditions  for  unique 
interpolation  of  polynomials 

Corollary  (2.4) 

Theorem  (2.8) 

rectangular 

necessary  conditions  for  unique 
reconstruction  of  BLP  signals 

Corollary  (2.5) 

Theorem  (2.7) 

rectangular 

recovery  of  bandlimited  signals 
of  a  given  order  from  an 
arbitrary  number  of  projections 

Table  2.1:  Summary  of  the  theorems  and  corollaries  of  Chapter  2. 


Chapter  3 


Reconstruction  Algorithms  for  the 
Semi-implicit  Sampling  Approach 

In  Chapter  2,  we  derived  a  spectrum  of  semi-implicit  sampling  strategies  for  a  variety  of 
multidimensional  reconstruction  problems.  Our  main  goal  in  this  chapter  is  to  derive  recon¬ 
struction  algorithms  for  these  sampling  schemes.  Although  most  of  the  examples  in  this  chapter 
are  related  to  the  problem  of  reconstruction  from  multiple  level  threshold  crossings,  the  majority 
of  our  proposed  algorithms  can  also  be  applied  to  the  other  reconstruction  problems  described 
in  section  (2.3). 

Our  first  approach  to  signal  reconstruction,  which  involves  solving  a  linear  system  of  equa¬ 
tions,  is  included  in  section  (3.1).  Although  this  approach  has  the  potential  of  being  numerically 
stable,  it  is  computationally  intensive  and  requires  large  amounts  of  storage.  Our  second  ap¬ 
proach  which  is  described  in  section  (3.2),  is  recursive  and  is  based  on  the  proofs  of  Theorems 
(2.3)  and  (2.4).  The  recursive  algorithm  requires  less  storage  and  computation.  However,  since 
it  computes  the  coefficients  recursively,  small  errors  in  the  initial  steps  of  the  algorithm  are 
propagated  and  magnified  in  the  final  steps  of  the  algorithm.  Thus,  as  we  expect,  unless  the 
interpolation  problem  at  hand  is  extremely  well  conditioned,  the  recursive  algorithm  does  not 
result  in  satisfactory  reconstruction.  In  section  (3.3),  we  will  derive  a  line  by  line  reconstruction 


algorithm  which  is  somewhat  similar  to  the  recursive  algorithm  of  section  (3.2),  but  is  free  of 
its  numerical  instabilities.  The  basic  idea  behind  line  by  line  reconstruction  is  to  recover  the 
one-dimensional  signal  associated  with  each  of  the  sampling  lines  and  then  to  interpolate  the 
original  signal  from  these  one-dimensional  functions.  We  will  propose  two  strategies  for  the 
former  part  of  the  algorithm.  The  first  strategy  is  non-iterative;  It  can  be  used  for  a  variety 
of  reconstruction  problems  utilizing  the  semi-implicit  sampling  strategy  The  second  strategy, 
which  is  iterative,  can  only  be  applied  to  the  problem  of  reconstruction  from  function  cross¬ 
ings.  In  addition,  the  number  of  interpolation  points  required  for  each  of  these  two  line  by  line 
reconstruction  strategies  exceeds  the  theoretical  requirements  of  Chapter  2.  The  main  features 
of  the  iterative  algorithm  are  its  robustness,  speed,  and  low  storage  requirements. 

For  all  the  reconstruction  examples  of  this  chapter,  the  location  of  crossings  are  specified 
to  16  digits.  The  quantization  issues  of  some  of  the  reconstruction  methods  described  in  this 
chapter  will  be  investigated  further  in  Chapter  5. 

3.1  Linear  Least-Squares  Approach 

The  most  straightforward  approach  to  the  reconstruction  of  signals  via  the  semi-implicit  sam¬ 
pling  strategy  is  to  solve  a  linear  system  of  equations.  Thus,  in  reconstructing  a  (21V  +  1)  x 
(21V  +  1)  BLP  signal  from  n(  thresholds,  we  must  first  find  a  set  of  sampling  lines  whose  in¬ 
tersections  with  the  threshold  contours  satisfy  the  theoretical  requirements  of  Chapter  2,  Then 

these  M  >  IV2  samples  are  used  to  solve  the  following  linear  system  of  equations: 

N  N 

Y  F[k1,k2)ei2ir{k'x>  +  =  U  1  <  t  <  n(,  1  <  j  <  M 

*,=-n  k2=-yv 

The  required  number  of  reconstruction  samples  depends  on  the  sampling  lines  and  is  at  least 
as  large  as  the  number  of  Fourier  coefficients.  As  we  will  see  in  section  (5.1.2),  increasing  the 


number  of  reconstruction  samples  improves  the  condition  number  of  the  reconstruction  problem 
at  hand. 

Examples  of  reconstruction  of  the  31  x  31  eye  picture  shown  in  figure  (2.10)  from  8  of 
its  level  crossings  are  shown  in  figures  (3.1)  and  (3.2).  Figure  (3.1)  shows  the  8  level  non- 
uniform  amplitude  quantized  version  of  figure  (2.10)  and  the  threshold  contours  associated 
with  these  levels.  Since  the  intensity  value  of  the  original  picture  lies  between  0  and  256,  the 
thresholds  were  also  chosen  in  this  range.  Figure  (3.2),  shows  the  reconstructed  versions  of 
the  original  eye  picture,  via  different  sets  of  sampling  lines.  More  specifically,  31  sampling 
lines  of  unit  slope  were  used  for  the  reconstruction  shown  in  figure  (3.2a),  11  sampling  lines 
of  slope  3  were  used  for  the  reconstruction  shown  in  figure  (3.2b),  and  15  sampling  lines  with 
slopes  1,2,  ±|,±3,±5,±4,  ±j,5  were  used  for  the  reconstruction  shown  in  figure  (3.2c).  In 
all  the  above  examples,  the  distributions  obtained  from  the  intersections  of  sampling  lines  and 
threshold  contours  satisfied  the  theoretical  requirements  of  Theorems  (2.3),  (2.4),  and  (2.3).  To 
improve  the  robustness,  we  used  all  the  intersections  of  sampling  lines  and  threshold  contours 
as  reconstruction  samples,  and  applied  QR  decomposition  [l]  to  solve  the  overdetermined  linear 
systems  of  equations  associated  with  each  example.  More  specifically,  the  number  of  samples 
used  for  reconstructions  shown  in  figures  (3.2a),  (3.2b),  and  (3.2c)  were  1475,  1266,  and  1548 
respectively. 

Robustness  of  the  linear  least-squares  approach  is  highly  dependent  upon  the  actual  al¬ 
gorithm  used  for  solving  the  overdetermined  system  of  equations.  Of  the  many  algorithms 
available,  we  chose  the  QR  decomposition  primarily  for  stability  reasons.  Unfortunately,  the 
most  stable  algorithms,  such  as  QR,  are  also  very  computation  intensive.  For  instance,  solving 
an  m  x  p  overdetermined  system  of  equations  requires  mp2  floating  point  operations  (flops). 


(a)  (b) 

Figure  3.1:  (a)  8  level  non  uniform  amplitude  quantized  ver¬ 
sion  of  the  original  eye  picture;  ( b )  Threshold  contours  asso¬ 
ciated  with  the  8  level  crossings. 

Thus,  assuming  that  the  nuuibe.  of  interpolation  points  used  is  of  the  same  order  of  magnitude 
as  the  number  of  unknown  Fourier  series  coefficients,  for  an  image  with  N  x  N  region  of  support 
in  the  Fourier  domain,  the  number  of  flops  required  by  the  QR  decomposition  algorithm  is  pro¬ 
portional  to  N 6.  The  other  major  drawback  of  the  linear  least-squares  approach  is  its  storage 
requirements.  The  storage  required  for  an  image  with  N  x  N  region  of  support  in  the  Fourier 
domain  is  N*.  Indeed,  storage  was  our  main  limiting  factor  in  terms  of  the  largest  image  we 
could  process  using  the  QR  decomposition. 

One  way  to  get  around  the  storage  problem  is  to  use  iterative  algorithms  such  as  the 
conjugate  gradient  method.  In  these  algorithms,  the  M  x  N2  matrix  associated  with  the  M 
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samples  of  a  signal  with  N  x  N  region  of  support  in  the  Fourier  domain  need  not  be  stored, 
since  in  each  step  of  the  iteration  we  only  need  to  multiply  it  implicitly  by  a  vector.  The 
conjugate  gradient  algorithm  needs  to  operate  on  symmetric  matrices.  Thus,  in  order  to  solve 
the  overdetermined  system  of  equations  given  by : 

A  x  =  b 

we  need  to  premultiply  both  sides  of  the  above  equation  by  AT  to  get: 

At  A  x  =  AT  b 

Since  the  condition  number  of  AT  A  is  the  square  of  the  condition  number  of  A,  and  the  conver¬ 
gence  rate  of  the  algorithm  is  proportional  to  the  condition  number  of  AT  A  [1],  if  our  original 
interpolation  problem  is  slightly  ill-conditioned  (i.e.,  A  has  rather  large  condition  number  ), 
then  it  will  take  many  iterations  before  the  algorithm  converges.  Thus,  the  conjugate  gradient 
approach  is  only  appropriate  for  situations  where  our  interpolation  problem  is  well-conditioned. 
Geometrically  speaking,  this  corresponds  to  the  situation  in  which  our  samples  are  more  or  less 
evenly  distributed  in  the  image.  We  have  found  experimentally  that  reconstruction  of  images 
from  multiple  level  threshold  crossings  via  the  conjugate  gradient  algorithm  is  unsuccessful, 
when  the  number  of  thresholds  is  small.  As  an  example,  consider  figure  (3.3b)  which  shows 
the  reconstructed  version  of  the  original  eye  picture  from  8  of  its  thresholds  crossings  on  lines 
of  slope  1  via  the  conjugate  gradient  algorithm.  The  intersections  of  threshold  crossings  and 
the  sampling  lines  ate  shown  in  figure  (3.3a).  As  is  apparent,  the  quality  of  reconstruction  in 
various  areas  of  the  image  is  highly  dependent  on  the  density  of  the  available  crossings  in  the 
region. 

Although  the  conjugate  gradient  algorithm  does  not  perform  satisfactorily  for  reconstruction 


(a) 


(b) 


Figure  3.3:  (a)  Intersections  of  lines  of  unit  slope  with  the  8 
hvel  crossings  shown  in  figure  (3.1b);  (b)  Reconstructed  ver¬ 
sion  of  the  eye  picture  via  the  conjugate  gradient  algorithm 
using  the  intersections  of  8  level  crossings  with  lines  of  unit 
slope. 


from  level  crossings,  it  can  be  used  successfully  for  reconstruction  problems  in  which  the  samples 
are  evenly  distributed  across  the  image.  Recall  from  section  (2.3.1)  that  the  sinusoid  crossings  of 
BLP  signals  lie  in  between  parallel  equidistant  lines  at  which  the  sinusoid  assumes  its  maximum 
and  minimum  values.  Thus,  we  would  expect  the  conjugate  gradient  algorithm  to  perform 
satisfactorily  for  the  problem  of  reconstruction  from  sinusoid  crossings1.  Figure  (3.4)  shows 
the  one  bit  amplitude  quantized  version  of  the  original  31  x  31  eye  picture  with  the  sinusoidal 
function 

'Note  that  even  distribution  of  the  samples  results  in  stable  reconstruction  of  the  signal  pi u*  timuoid,  and 
not  necessarily  the  signal  itself.  For  instance,  if  the  dynamic  range  of  the  sinusoid  is  much  larger  than  that  of 
tha  signal,  we  would  expect  the  quality  of  reconstruction  of  the  signal  by  itself  to  be  rather  poor. 
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128  {1  +  cos(2x(16x  +  16y)j} 

Figure  (3.5)  shows  the  reconstructed  version  of  the  eye  picture  via  the  conjugate  gradient 
algorithm  with  20  iterations.  The  reconstruction  samples  were  chosen  to  be  a  subset  of  all  the 
intersections  of  the  31  lines  of  unit  slope  with  the  crossings,  satisfying  the  minimum  distribution 
requirements  of  Corollary  (2.1). 

Our  second  example  involves  reconstruction  of  the  63  x  63  picture  shown  in  figure  (3.6) 
which  will  be  referred  to  as  cman.  One  bit  amplitude  quantized  version  of  the  sum  of  original 
cman  picture  and  the  sinusoidal  function: 

128{1  +  cos[2jt(32x  +  32y)]> 

is  shown  in  figure  (3.7).  Its  reconstructed  versions  via  the  conjugate  gradient  algorithm  with  20 
iterations  are  shown  in  figure  (3.8).  Figure  (3.8a)  was  reconstructed  from  ail  the  intersections  of 
63  equidistant  lines  of  slope  one  with  sinusoid  crossings,  whereas  a  smaller  subset  of  the  inter¬ 
sections  satisfying  the  theoretical  requirements  of  Theorem  (2.3)  were  chosen  for  reconstruction 
of  figure  (3.8b).  Thus,  all  the  sampling  lines  in  the  former  case  have  64  points,  whereas  in  the 
latter  case  there  are  sampling  lines  with  as  few  points  as  1  or  3.  Indeed,  the  artifacts  shown 
in  figure  (3.8b)  are  due  to  the  fact  that  some  sampling  lines  contain  many  fewer  interpolation 
points  than  other  lines.  To  summarize  this  section,  we  have  proposed  the  application  of  two 
major  least-squares  algorithms  for  solving  the  linear  system  of  equations  associated  with  our 
reconstruction  problem.  The  QR  decomposition  approach  is  extremely  robust  and  can  be  used 
for  a  variety  of  reconstruction  problems  including  the  ones  in  which  the  samples  are  not  evenly 
distributed  across  the  image.  The  drawback  of  the  QR  decomposition  approach  is  that  it  is 
computationally  intensive  and  its  storage  requirements  are  rather  large.  The  conjugate  gra- 
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Figure  3.5:  Reconstruction  of  the  eye  picture  from  its  sinusoid  crossings  via  the  conjugate 
gradient  algorithm.  Reconstruction  samples  were  chosen  from  a  subset  of  all  the  intersections 
of  lines  of  unit  slope  with  the  crossings,  satisfying  the  minimum  distribution  requirements  of 
Theorem  (3.6). 
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dient  algorithm  needs  less  storage  and  converges  rather  quickly  for  well  conditioned  problems 
such  as  reconstruction  from  sinusoid  crossings.  However,  its  slow  convergence  properties  for 
ill  conditioned  problems  limits  ts  applicability  to  reconstruction  from  multiple  level  threshold 
crossings. 


3.2  Recursive  Approach 


As  mentioned  in  Chapter  2,  although  the  algebraic  geometric  results  of  section  (2.2.3)  are 
considerably  more  general  than  those  of  sections  (2.2.1)  and  (2.2.2),  the  inherent  advantage  of 
the  results  in  the  latter  sections  lies  in  the  existence  of  recursive  algorithms  to  carry  out  the 
reconstruction  proposed  by  their  theoretical  results2 .  These  algorithms  are  essentially  based  on 
the  constructive  proofs  of  Theorems  (2.3)  and  (2.4)  which  provide  us  with  sufficient  conditions 
for  unique  recovery  of  BLP  signals  from  their  samples  on  lines  with  fixed  positive  integer  slope. 
These  proofs  are  constructive  in  the  sense  that  unique  recovery  is  shown  by  actually  deriving 
a  reconstruction  strategy.  Since  detailed  description  of  the  recursive  algorithms  are  included 
in  Chapter  2  and  Appendix  (A),  we  will  briefly  describe  the  recursive  algorithm  based  on  the 
proof  of  Theorem  (2.3)  and  show  a  few  examples  of  reconstruction. 

Theorem  (2.3)  of  section  (2.2.1)  provides  us  with  the  distribution  of  points  on  lines  passing 
through  the  origin,  which  results  in  unique  recovery  of  polynomials 


p{w,z)  =  YLY,  a(‘-j)  w>  z 3  (31) 

«=0  ;= 0 

More  specifically,  it  states  that  for  unique  reconstruction  of  a  polynomial  of  the  above  form,  we 

sIt  might  also  be  possible  to  derive  a  recursive  proof  for  the  results  of  section  (2.2.3).  However,  since  such  a 
proof  would  be  considerably  more  complex  than  the  algebraic  geometric  one,  and  its  resulting  recursive  algorithm 
is  likely  to  be  rather  unstable  (for  reasons  which  we  will  discuss  in  this  section),  we  have  decided  not  to  presue 
the  recursive  approach  for  the  theoretical  results  of  section  (2.2.3). 
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with  2t  4-  1  points  on  the  tth  line.  Defining  the  variables  b[k^  associated  with  the  fcth  sampling 

line  in  terms  of  its  slopes  a*  and  the  polynomial  coefficients  a(t,j)  in  the  following  way: 

'  « 

Y  a('  ~  m>  >  0<i<n 

6t(*J  =  m=n°  (3.2) 

Y  a(*  “  m>  m)ajT  >  n  <  i  <  2n 

m—t—n 

we  conclude  that  samples  on  the  fcth  line  satisfy: 

Y  bik)w'  -  p{w,aku>)  -  Ybtk)w'  ~  Y.  b\k)w'  (3-3) 

i=j+1  «=0  i=2n-4 

The  recursive  algorithm  essentially  consists  of  2n  +  1  steps  corresponding  to  s  =  -1,  ...,2n-  1 
in  the  above  equation.  More  specifically,  at  the  sth  stage,  the  algorithm  uses  the  computed 
values  of  and  bi^  for  0  <  k  <  n,  together  with  the  2(n  -  s)  -  1  samples  on  the  line 

ln-3-i ,  in  order  to  find  b\n  3  for  s  +  1  <  t  <  2n  -  s  -  1  by  solving  a  linear  system  of 
equations  indicated  by  equation  (3.3).  These  values  are  then  used  in  equation  (3.2)  to  find  the 
coefficients  {a(i,j)  |  i  +  j  =  s  +  1,  2n  -  s  -  1}  by  solving  a  second  set  of  linear  equations. 
Finally,  the  coefficients  {a(t,j)  |  t  +  j  =  s  +  l,2n  -  s  —  1}  are  used  to  find  an<^  bi+ 1 

for0<fc<n-s-2  via  equation  (3.2)  for  the  next  step  of  the  algorithm. 

From  the  above  description,  we  see  that  the  polynomial  coefficients  are  found  recursively  and 
that  the  coefficients  found  in  initial  steps  are  used  to  compute  the  ones  in  the  later  stages.  Thus, 
as  we  might  expect,  small  errors  in  finding  the  coefficients  in  the  initial  steps  are  propagated 
and  magnified  in  the  final  steps  of  the  algorithm.  Indeed,  this  turns  out  to  be  a  major  source 
of  instability  for  tl.e  recursive  algorithm. 


Since  BLP  signals  can  be  expressed  as  polynomials,  the  above  recursive  algorithm  can  be 
used  for  reconstructing  signals  from  their  samples  on  lines  of  unit  slope.  In  applying  this  al¬ 
gorithm  to  the  problem  of  reconstruction  from  multiple  level  crossings,  we  have  found  that 
unless  the  interpolation  points  are  more  or  less  evenly  distributed  across  the  image,  the  al¬ 
gorithm  fails  to  recover  the  signal  successfully.  For  instance,  the  algorithm  is  unsuccessful  in 
reconstructing  the  31  x  31  eye  picture  of  figure  (2.10)  from  its  8  level  crossings  shown  in  figure 
(3.1b),  even  though  it  is  capable  of  reconstructing  it  from  its  sinusoid  crossings  shown  in  figure 
(3.4).  This  is  because,  as  we  mentioned  earlier,  the  sinusoid  crossings  of  the  signal  are  almost 
uniformly  distributed  across  the  image.  The  reconstructed  version  of  the  eye  picture  from  its 
crossings  with  128{1  +  cos[ 2jt  (16x  +  16y)]}  is  shown  in  figure  (3.9).  The  reconstruction 
samples  were  chosen  to  be  a  subset  of  all  the  intersections  of  31  equidistant  lines  of  unit  slope 
with  the  crossings,  satisfying  the  minimum  required  distribution  of  Theorem  (2.3). 

The  reconstructed  version  of  the  cman  picture  from  its  sinusoid  crossings  shown  in  figure 
(3.7),  is  shown  in  figure  (3.10).  The  reconstruction  samples  were  chosen  to  be  a  subset  of 
the  intersections  of  63  equidistant  lines  of  unit  slope  and  the  sinusoid  crossings,  satisfying  the 
distribution  requirements  of  Theorem  (2.3).  Notice  that  unlike  reconstruction  via  conjugate 
gradient  algorithm  shown  in  figure  (3.8b),  there  are  no  artifacts  in  figure  (3.10)  due  to  unequal 
number  of  samples  on  different  sampling  lines.  This  indicates  that  the  recursive  algorithms  are 
somewhat  more  stable  than  the  conjugate  gradient  algorithm.  In  fact,  the  recursive  algorithm 
described  in  this  section  is  robust  enough  for  reconstruction  from  sinusoid  crossings  of  larger 
images  such  as  the  one  shown  in  figure  (3.11).  Figure  (3.11)  which  will  be  referred  to  as 
vegas,  has  127  x  127  Fourier  coefficients.  Its  crossings  with  128  (1  +  cos [2x  (64x  +  64y)]}  are 
shown  in  figure  (3.12)  and  its  reconstructed  version  is  shown  in  figure  (3.13).  The  reconstruction 
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Finally,  let  us  briefly  discuss  the  efficiency  and  storage  requirements  of  the  recursive  algo¬ 
rithm,  described  in  this  section.  As  we  mentioned  earlier,  to  find  the  coefficients  of  a  polynomial 
given  by  equation  (3.1),  the  recursive  algorithm  goes  through  2n  +  1  steps  and  in  each  step  it 
has  to  solve  two  Van  der  Monde  systems  of  linear  equations.  More  specifically,  in  the  kth  step, 
the  size  of  these  systems  are  (2n  —  k)  x  (2n  —  fc)  and  (k  +  2)x  (k  +  2).  Although  there  are  ways  of 
solving  a  p  x  p  Van  der  Monde  system  of  equations  with  0(p2)  flops,  these  methods  are  usually 
unstable.  Stable  techniques,  on  the  other  hand,  require  0(p3)  floating  point  operations.  Thus 
the  total  number  of  operations  required  for  reconstruction  of  a  signal  with  N  x  N  region  of 
support  in  the  Fourier  domain  is  of  the  order  of  0(N*).  Furthermore,  the  storage  requirements 
of  the  recursive  algorithm  is  of  the  order  of  N2. 

Comparing  characteristics  of  the  recursive  algorithm  with  the  QR  decomposition  we  find 
that  the  latter  requires  more  computation  and  space  but  is  also  more  stable.  However,  as  we  saw 
in  figures  (3.8b)  and  (3.10),  the  recursive  approach  is  more  stable  than  the  iterative  approach 
of  the  conjugate  gradient  algorithm.  Furthermore,  the  recursive  algorithm  can  only  be  used 
for  reconstruction  from  samples  on  lines  of  integer  slope  with  all  the  lines  having  equal  slopes, 
whereas  the  linear  least-squares  approaches  of  section  (3.1)  are  capable  of  reconstructing  from 
samples  on  lines  with  different  positive  or  negative  rational  slopes.  However,  the  least-squares 
and  recursive  methods  are  similar,  in  the  sense  that  they  can  both  be  used  for  the  more  general 
problem  of  reconstruction  from  non-uniformly  distributed  samples,  as  opposed  to  the  specific 
problem  of  reconstruction  from  level  crossings.  As  we  will  see  in  the  next  section,  this  is  not 
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Figure  3.9:  Reconstructed  version  of  the  eye  picture  via  the  recursive  algorithm  from  its  sinusoid 
crossings.  The  reconstruction  samples  were  chosen  to  be  a  subset  of  ail  the  intersections  of  31 
equidistant  lines  of  unit  slope  with  the  crossings,  satisfying  the  minimum  required  distribution 
of  Theorem  (3.6). 


Figure  3.10:  Reconstructed  version  of  the  cman  picture  from  its  sinusoid  crossings  via  the 
recursive  algorithm.  The  reconstruction  samples  were  chosen  to  be  a  subset  of  all  the  intersec¬ 
tions  of  63  equidistant  lines  of  unit  slope  with  the  crossings,  satisfying  the  minimum  required 
distribution  of  Theorem  (3.6). 


Figure  3.11:  The  original  vegas  picture  with  127  x  127  region 
of  support  in  the  Fourier  domain. 


the  case  for  the  iterative  version  of  our  line-by-line  reconstruction  algorithm. 


3.3  Line  by  Line  Reconstruction  Approach 


As  we  saw  in  the  previous  section,  the  major  drawback  of  the  recursive  approach  to  recon¬ 


struction  is  its  numerical  instability,  which  is  a  direct  consequence  of  computing  the  Fourier 


coefficients  recursively.  More  specifically,  the  recursive  algorithm  finds  the  coefficients  in  the 


kth  step  by  using  the  value  of  the  coefficients  it  has  computed  in  the  first  (k  -  1)  steps.  Thus,  as 


we  might  expect,  small  errors  in  computing  the  coefficients  in  the  initial  steps  are  propagated 


and  magnified  in  the  final  steps  of  the  algorithm. 


Our  main  goal  in  this  section  is  to  propose  a  new  algorithm,  similar  in  style  to  the  recursive 
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Figure  3.12:  One  bit  amplitude  quantized  version  of  the  sum  of  vegas  picture  and 
128{1  +  cos[ 2x  (64z  +  64y)]}. 


Figure  3.13:  Reconstructed  version  of  the  vegas  picture  from  its  sinusoid  crossings  via  the 
recursive  algorithm.  The  reconstruction  samples  were  chosen  to  be  a  subset  of  all  the  intersec¬ 
tions  of  127  equidistant  lines  of  unit  slope  with  the  crossings,  satisfying  the  minimum  required 
distribution  of  Theorem  (3.6). 
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approach,  which  is  free  of  its  numerical  instabilities.  The  line  by  line  reconstruction  strategy  we 
are  going  to  describe  in  this  section  exploits  the  fact  that  the  one-dimensional  signal  obtained 
by  sampling  our  two-dimensional  signal  along  a  line  with  rational  slope  is  bandlimited  and 
periodic.  More  specifically,  consider  a  two-dimensional,  BLP  signal  with  (2 N  +  1)  x  (2 N  +  1) 
region  of  support  in  the  Fourier  domain  given  by 

f{*,v)  =  E  E  F{kukt)  e**”k'  <?'*** 

ki=-N  k2=-N 


If  it  is  sampled  along  the  line: 


my  =  nx  +  0 


the  one-dimensional  signal  is  periodic  and  has  2iV(m-fn)  +  l  Fourier  series  coefficients.  Thus 
2 N(m  +  n)  +  1  samples  of  this  one-dimensional  signal  will  enable  us  to  uniquely  specify  it3. 
Since  our  semi-implicit  line  sampling  strategy  of  Chapter  2  consists  of  using  points  on  lines  of 
various  rational  slopes,  an  obvious  way  to  recover  a  two-dimensional  signal  would  be  to 

1.  Reconstruct  each  of  the  one-dimensional  signals  corresponding  to  the  sampling  lines  sep¬ 
arately. 

2.  Recover  the  two-dimensional  signal  using  the  one-dimensional  ones. 

We  now  need  to  address  a  few  issues  with  regard  to  the  above  reconstruction  strategy.  The 
first  issue  is  whether  or  not  the  sampling  distribution  proposed  by  the  theoretical  results  of 
Chapter  2  provide  us  with  enough  interpolation  points  to  recover  the  one-dimensional  signal 
associated  with  each  of  the  sampling  lines.  As  it  turns  out,  the  sampling  requirements  of 

the  above  line  by  line  reconstruction  scheme  are  more  stringent  than  those  of  Corollary  (2.3), 

3The  fact  that  the  one-dimensional  signal  associated  with  a  sampling  line  with  rational  slope  is  bandlimited 
and  periodic  was  originally  shown  by  Merseraeu  [2j. 
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which  is  the  most  general  result  of  Chapter  2.  For  instance,  in  situations  where  all  of  the 
sampling  lines  have  integer  slopes  p,  the  theoretical  results  of  Chapter  2  require  2(p+  1 )  jV  -  2 pi 
interpolation  points  on  the  ith  line,  whereas  the  above  line  by  line  reconstruction  strategy  needs 
2 (p+  1  )N  samples  on  all  the  lines.  In  other  words,  for  line  by  line  reconstruction,  the  number 
of  samples  required  on  each  of  the  lines  is  equal  to  the  maximum  number  of  points  required  by 
Theorem  (2.4)  on  any  one  of  them.  Since  the  line  by  line  approach  requires  more  samples  than 
the  minimum  specified  by  the  theoretical  results  of  Chapter  2,  we  would  expect  that,  for  the 
problem  of  reconstruction  from  level  crossings,  more  thresholds  are  needed  to  generate  these 
samples.  Indeed,  as  we  will  see  in  the  examples  of  this  section,  while  the  eye  picture  of  figure 
(2.10)  can  be  reconstructed  from  samples  of  its  8  level  crossings  via  the  linear  least-squares 
approach,  14  level  crossings  are  needed  for  its  line  by  line  reconstruction. 

The  second  issue  which  needs  to  be  addressed  has  to  do  with  the  recovery  of  the  two- 
dimensional  signal  under  consideration,  once  all  the  one-dimensional  signals  associated  with 
the  sampling  lines  have  been  reconstructed.  In  situations  where  all  the  nodes  of  the  (21V  + 1)  x 
(2 N  +  1)  square  grid  of  the  two-dimensional  signal  are  located  on  the  sampling  lines,  the  value 
of  the  signal  on  the  grid  can  be  computed  easily  from  the  one-dimensional  signals.  Equally 
spaced  sampling  lines  of  unit  slope  is  an  example  of  such  a  situation.  For  the  more  general  case, 
we  can  proceed  as  follows: 

1.  Use  the  one-dimensional  signals  to  find  their  intersections  with  (2 N  +  1)  equally  spaced 
horizontal  or  vertical  lines.  For  the  most  general  sampling  scenario,  described  by  Corollary 
(2.3),  we  can  show  that  each  horizontal  or  vertical  line  intersects  our  sampling  lines  in  at 
least  (2 N  +  1)  points.  More  specifically,  if  inequality  (2.40)  of  Corollary  (2.3)  is  satisfied, 
then  each  of  the  (21V  +  1)  horizontal  lines  intersect  our  sampling  lines  at  (21V  +  1)  points; 


On  the  other  hand,  if  inequality  (2.41)  is  satisfied,  then  each  of  the  (2N  +  1)  vertical  lines 
intersect  our  sampling  lines  at  ( 2N  +  1)  points. 

2.  Since  the  one-dimensional  signal  associated  with  each  horizontal  or  vertical  line  is  a  BLP 
signal  with  (2 N  +  1)  Fourier  coefficients,  they  are  reconstructible  from  their  (2 N  +  1) 
samples,  obtained  in  the  previous  step. 

3.  For  each  of  the  horizontal  or  vertical  lines,  find  the  values  of  the  signal  on  (2 N  + 1)  equally 
spaced  points.  These  (2 N  +  1)  points  on  each  of  the  (2 N  -I-  1)  horizontal  or  vertical  lines 
correspond  to  samples  of  the  two-dimensional  signal  on  a  (2 N  +  1)  x  (2N  +  1)  square  grid. 

The  third  issue  we  need  to  address  is  the  recovery  of  the  one-dimensional  signals  associated 
with  the  sampling  lines.  As  we  mentioned  earlier,  the  sampling  requirements  of  line  by  line 
reconstruction  could  be  more  stringent  than  those  required  by  the  theoretical  results  of  Chapter 
2.  Since  the  one-dimensional  signal  along  the  ith  sampling  line: 

l,  :  m,  y  =  tij  i  +  ft 

is  bandlimited  and  periodic  with  2 N(nii  -f-  n,)  1  Fourier  harmonics,  theoretically  speaking, 

2N(mi  +  rij)  +  1  samples  of  it  are  sufficient  to  reconstruct  it  uniquely.  We  can  find  its  Fourier 
series  coefficients  by  solving  a  (2 N(rrii  +  rq)  +  1)  x  (2 N(mi  +  rq)  +  1)  linear  system  of  equations, 
or  use  Lagrange  interpolation  to  compute  its  samples.  An  example  of  reconstruction  of  the 
original  eye  picture  shown  in  figure.  (2.10)  from  intersections  of  14  of  its  threshold  crossings 
with  lines  of  unit  slope  is  shown  in  figure  (3.14).  In  this  example,  the  one-dimensional  signals 
associated  with  the  sampling  lines  were  found  by  solving  linear  systems  of  equations. 

For  the  specific  problem  of  reconstruction  from  level  crossings,  we  can  propose  an  iterative 
algorithm  which  uses  all  the  intersections  of  the  level  crossings  with  sampling  lines  to  reconstruct 
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Figure  3.14:  Reconstruction  of  the  original  eye  picture  from 
intersections  of  lines  of  unit  slope  with  14  of  its  level  crossings. 
Non-iterative  line  by  line  reconstruction  was  used. 


the  one-dimensional  signals  associated  with  the  sampling  lines.  We  will  devote  the  remaining 


part  of  this  section  to  this  iterative  algorithm. 


3.3.1  Iterative  Approach  to  Line  by  Line  reconstruction 


In  this  section,  we  will  describe  an  iterative  algorithm  for  the  reconstruction  of  one-dimensional 


signals,  associated  with  sampling  lines,  from  their  intersections  with  level  crossings.  As  we 


will  see,  similar  to  the  least-squares  and  recursive  approaches  of  sections  (3.1)  and  (3.2),  the 


iterative  approach  can  be  applied  both  to  the  specific  problem  of  reconstruction  from  function 


crossings  and  to  the  more  general  problem  of  reconstruction  from  non-uniformly  distributed 


samples. 
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The  basic  idea  behind  iterative  reconstruction  is  to  iterate  between  space  and  frequency 
domains  by  imposing  appropriate  constraints  on  the  signal  in  each  domain.  The  frequency 
domain  constraint  is  due  to  the  fact  that  the  one-dimensional  signal,  obtained  by  sampling 
our  two-dimensional  signal  along  a  line  with  rational  slope,  is  bandlimited  and  periodic.  The 
space  domain  constraint  can  be  derived  using  the  level  crossing  information.  More  specifically, 
if  the  locations  of  all  the  crossings  of  a  one-dimensional,  continuous,  BLP  signal  g(z)  with  p 
thresholds  <  £2  •••  <  tp  are  known,  and  the  signal  is  known  to  be  in  the  range  \to,tp+i], 
then  for  any  arbitrary  zq,  we  can  deduce  the  range  into  which  j(zo)  falls.  That  is,  given  zo,  we 
can  find  1  such  that 

t,  <  g(za)  <  t,+l 

The  above  process  is  shown  pictorially  in  figure  (3.15).  In  order  to  derive  the  range  information 
for  any  point  on  a  one-dimensional  signal,  we  need  a  minimum  of  two  samples  corresponding 
to  two  different  thresholds.  This  has  to  do  with  the  fact  that  the  samples,  corresponding  to 
one  threshold,  are  incapable  of  resolving  the  sign  ambiguity  of  the  one-dimensional  signal  under 
consideration. 

Using  the  above  ideas,  we  propose  the  following  algorithm  for  reconstructing  a  one-dimensional, 
continuous,  BLP  signal,  j(z)  with  period  1  and  2 N  +  1  Fourier  harmonics  from  all  of  its  level 
crossings  with  thresholds 

1.  Deduce  the  range  of  intensity  for  M  >  2 N  +  1  equally  spaced  points 

1  2  M  -  1 

’  M'  M’  ’  M 

That  is,  for  each  of  the  above  points,  find  the  threshold  t,  such  that  the  actual  value  of  g 
at  that  point  is  larger  than  U  and  smaller  than  t,+1.  Thus,  if  ti(n)  and  tu(n)  denote  the 


•  points  with  intensity  <  t, 

©  points  with  intensity  >  ti  and  <  t 
0  points  with  intensity  >  t 


O  points  with  intensity  >  t 
•  points  with  intensity  <  t 


Figure  3.15:  (a)  Deriving  the  space  domain  constraint  for  recovery  of  a  one-dimensional  sig¬ 
nal  via  the  iterative  algorithm;  (b)  Ambiguous  range  information  via  samples  of  the  contour 
associated  with  one  threshold. 
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lower  and  upper  bound  for  the  nth  point,  we  have: 


ti[n)  <  g(-)  <tu(n)  ,  n  =  0,  (3.4) 

2.  Let  £^(^7)  denote  the  value  of  g  in  the  fth  iteration  at  the  point  ^  .  Then,  the  initial 


guess  becomes: 


»<»»(£>  =  .  „  -  0 . M-  1 


As  we  will  see  later,  the  convergence  of  the  algorithm  is  independent  of  the  initial  guess. 


3.  Take  the  discrete  Fourier  transform  (DFT)  of  g  in  the  /th  step: 


G(I>(fc)  =  DFT\gM{jj)\  ,  0  <  n,k  <  M 


4.  Impose  the  bandlimited  constraint: 


G(<+1)(fc)  = 


G(,)(fc)  ,  0<k<N,M-N<k<M 


,  elsewhere 


5.  Take  the  inverse  discrete  Fourier  transform  ( DFT  *)  of  G^l+l\k): 


=  DFT_1jG(i+,)(fc)]  ,  0  <  k,n  <  M 


6.  Impose  the  space  domain  constraint: 

9{l+l}(Zf)  ,  U(n)  <  t«{n) 

g[l+l) =  '  */(»)  ,  ?(1+,)(^)  <  U(n)  (3-5) 

,  5(i+1'(j&)  >  M") 

7.  If  all  the  M  points  of  <7^+IKs?)  satisfy  the  space  domain  constraint  to  within  some 
acceptable  error,  then  we  are  done.  Otherwise  go  to  step  (3)  and  repeat  steps  (3)  through 


«iLSl. 
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As  far  as  convergence  of  the  above  algorithm  goes,  we  can  show  that  the  mean  square  error 
between  the  unknown  signal,  which  satisfies  the  space  and  frequency  domain  constraints,  and 
the  solution  obtained  via  successive  iterations  decreases.  Furthermore,  as  we  will  see  in  the 
next  section,  we  can  apply  the  theory  of  projection  onto  convex  sets  (POCS)  to  prove  that  the 
algorithm  is  guaranteed  to  converge  to  a  solution,  which  satisfies  both  the  space  and  frequency 
domain  constraint,  if  such  a  solution  exists. 

3.3.2  Application  of  POCS  to  the  Iterative  Algorithm 

We  could  use  the  theory  of  projection  onto  convex  sets  to  establish  the  convergence  of  our 
algorithm.  As  we  will  see,  our  algorithm  is  a  special  case  of  a  more  general  iterative  scheme 
discussed  extensively  in  the  literature  [46,47] .  By  generalizing  our  algorithm,  we  can  not  only 
establish  its  convergence,  but  also  find  ways  of  improving  its  convergence  properties.  The  basic 
theorem  describing  the  general  iterative  method  can  be  stated  as  follows: 
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Theorem  3.1  (Youla  and  Webb  [46])  Let  H  be  a  Hilbert  space  with  elements  f,g,...,etc, 
a  zero  vector  </>  and  an  inner  product  (x,  y).  Furthermore,  let  Cq,  the  intersection  of  closed, 
convex  subsets  C j,  ...,Cm  of  H  be  given  by: 

Co  =  n£,C. 

be  non-empty.  Consider  the  composition  operator 

T  =  TmTm-i...T, 


where 

Ti  =  1  +  A i(Pt  -  1) 

and  Pt  is  a  projection  operator  onto  C,  .  Then,  for  every  x  €  H  and  every  choice  of  relaxation 
constants  Ai,...,Am  in  the  interval: 

0  <  A,  <  2  ,  1  <  i  <  m 

the  sequence  {Tnx}  converges  weakly  to  a  point  in  Cq.  If  H  is  finite  dimensional,  the  conver¬ 
gence  is  strong.  In  addition,  if  one  of  the  Ci ’s  say  Cm,  is  finite  dimensional,  by  setting 


•V 


A 
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we  are  guaranteed  of  strong  convergence  even  though  H  might  be  infinite  dimensional. 


Before  we  proceed  to  apply  the  above  result  to  our  reconstruction  problem,  let  us  briefly  go 
over  a  few  definitions.  The  sequence  {/*}  is  said  to  converge  strongly  to  /  if 


lim  ||  ft  ~  f  II  =  0 
*-—♦00 


The  sequence  {/*}  is  said  to  converge  weakly  to  /  if 


lim  {fk,g)  =  ( f,g )  Vy  <=  H 

fc— *oo 


Finally,  the  projection 


h  =  P,  f 


of  /  onto  C,  is  the  element  which  satisfies 


min  ||  /  -  y  ||2  = 

V6  C, 


/  ~  >l||J 


Our  iterative  algorithm  is  essentially  a  special  case  of  Theorem  (3.1)  where  the  Hilbert  space 
H  is  the  finite  dimensional  space  of  all  M  point  real  sequences  and  the  inner  product  between 
two  M  point  sequences  z(n)  and  y(n)  is  given  by. 

M-i 

(*.»)  =  Yi  I(n)v(n) 

n=0 

The  convex,  closed  subsets  of  H,  are  given  by  C\  and  C2.  C\  is  the  set  of  all  M  point  real, 
bandlimited  sequences  whose  DFT  is  0  in  the  range  N  <  k  <  M  -  N ,  and  C2  is  the  set  of 
all  M  point  real  sequences,  y(n),  which  satisfy  the  space  domain  constraint  derived  from  the 
level  crossing  information  as  described  in  the  first  step  of  our  algorithm.  That  is: 


*i(n)  ^  y(n )  <  tu(n)  0  <  n  <  M 


Pi 
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Note  that  for  C2  to  be  closed,  it  is  necessary  to  define  it  exactly  the  way  it  is  shown  above.  For 
instance,  if  in  inequality  (3.6),  the  <  were  replaced  with  <,  then  Ci  would  not  be  closed.  Pi 
is  the  bandlimiting  operator  which  projects  onto  the  set  C\ ,  and  P%  is  the  projection  operator 
which  imposes  the  space  domain  constraint  and  projects  onto  C 2.  Thus,  our  iterative  algorithm 
is  nothing  but  a  successive  application  of  the  operator  T  —  Pi  P2  to  an  arbitrary  initial 
condition,  which  we  chose  in  the  second  step  of  the  algorithm.  Therefore,  taking  into  account 
the  fact  that  H  is  a  finite  dimensional  vector  space,  and  applying  Theorem  (3.1)  with 

Ai  =  A2  =  1 

we  can  conclude  that  our  iterative  algorithm  converges  strongly  to  an  element  in  the  set  Co, 
which  in  this  case,  is  the  set  of  M  point  real  sequences  which  satisfy  both  the  space  and 
frequency  domain  constraints.  It  is  worthwhile  to  mention  that,  since  M  is  finite,  there  might 
be  many  M  point  real  sequences  in  the  set  Co.  However,  in  the  limit  as  M  — ►  00,  the  set  C0 
will  contain  exactly  one  sequence,  as  long  as  the  number  of  samples  on  a  given  line  exceeds  the 
number  of  Fourier  coefficients  of  the  one-dimensional  signal  associated  with  the  line  (because  of 
the  uniqueness  of  Fourier  series).  Under  these  circumstances,  H  will  not  be  finite  dimensional; 
thus  one  might  begin  to  worry  whether  strong  convergence  is  guaranteed  any  longer.  However, 
as  it  is  mentioned  in  Theorem  (3.1),  since  at  least  one  of  our  convex  closed  suDsets,  namely  Cj, 
is  finite  dimensional,  we  are  still  guaranteed  of  strong  convergence  as  long  as  we  set 

Ai  =  1 

Theorem  (3.1)  not  only  outlines  the  conditions  under  which  our  iterative  algorithm,  de¬ 
scribed  earlier  in  this  section,  converges,  but  also  generalizes  our  algorithm  by  introducing  the 
relaxation  parameters  Aj  and  A2.  We  can  control  the  convergence  properties  of  our  algorithm 


by  adjusting  the  values  of  the  A’s.  We  speak  of  under-relaxation  or  over-relaxation,  depending 
on  whether  0<A<lorl<A<2  respectively.  With  A  =  2,  we  have  another  special 
case,  referred  to  by  Motzkin  and  Schoenberg  [48]  els  the  reflection  method.  Values  of  A  outside 
the  interval  (0, 2]  disrupt  convergence. 

Convergence  of  the  generalized  iterative  algorithm  depends  on  the  “angle”  formed  at  the 
intersection  of  the  projection  sets,  C\  and  Ci\  a  wide  angle  between  the  sets  results  in  a  faster 
convergence  than  situations  where  the  sets  intersect  at  a  smaller  angle.  In  general,  one  can 
expect  that  the  smaller  the  intersection  of  the  sets  is,  the  smaller  the  angle  at  which  they 
intersect  will  be,  resulting  in  slower  convergence.  There  are  many  ways  to  accelerate  the 
convergence  of  this  type  of  iterative  algorithm.  The  simplest  such  technique  is  over- relaxation, 
which  essentially  consists  of  choosing  the  relaxation  parameter  between  1  and  2;  this  has  the 
effect  of  “opening  up”  a  small  angle  of  intersection  between  sets.  As  we  will  see  in  the  next 
section,  for  the  particular  problem  of  reconstruction  from  level  crossings,  convergence  rates  do 
not  present  serious  enough  problems  for  us  to  seek  other  means  of  accelerating  convergence. 

3.3.3  Examples  of  the  Iterative  Algorithm 
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An  example  of  reconstruction  of  the  original  31  x  31  eye  picture  of  figure  (2.10)  via  the 
iterative  algorithm  is  shown  in  figure  (3.16).  Intersections  of  14  level  crossings  with  32  lines  of 
unit  slope  were  used  to  recover  the  one-dimensional  signal  associated  with  the  sampling  lines. 
Although  our  theoretical  results  of  Chapter  2  only  require  31  sampling  lines,  to  be  able  to 
use  the  FFT  with  a  power  of  2  for  interpolation  from  one-dimensional  signals  to  the  square 
grid,  we  chose  32  equally  spaced  sampling  lines.  Thus  samples  of  the  two-dimensional  signal 
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Figure  3.16:  Reconstructed  version  of  the  original  eye  picture 
from  intersections  of  lines  of  unit  slope  with  14  level  crossings 
via  the  iterative  algorithm. 


on  a  32  x  32  grid  coincided  with  equally  spaced  points  on  the  sampling  lines.  The  relaxation 
parameters  were  chosen  to  be  A x  =  A  2  =  1,  and  the  number  of  equally  spaced  points  on  the 
sampling  lines,  M,  was  chosen  to  be  64. 

The  main  reason  for  choosing  such  a  large  number  of  thresholds  in  the  above  example  is  that 
the  number  of  crossings  on  a  given  line  must  be  larger  than  the  number  of  Fourier  harmonics 
of  the  one-dimensional  signal  associated  with  that  line.  However,  even  if  this  condition  is 
satisfied,  there  might  still  be  many  one-dimensional,  BLP  signals  satisfying  the  space  and 
frequency  domain  constraints,  as  long  as  M  is  finite.  Thus,  unique  reconstruction,  via  the 
iterative  algorithm,  is  possible  only  in  the  limit  as  M  — *  00,  if  the  number  of  crossings  on  each 
line  is  greater  than  or  equal  to  the  number  of  Fourier  series  coefficients  of  the  one-dimensional 


signal  associated  with  the  line.  Since  in  practice,  M  can  not  be  infinite,  the  solution  obtained  by 


the  iterative  algorithm  is  only  an  approximate  one,  regardless  of  the  number  of  crossings  used 
for  reconstruction.  Therefore,  the  question  which  arises  naturally  is  the  way  reconstruction 
quality  is  degraded  as  the  number  of  crossings  is  decreased  below  the  number  of  Fourier  series 
coefficients  of  the  one-dimensional  signal  to  be  recovered.  We  will  answer  this  question  more 
quantitatively  in  Chapter  5.  In  this  section,  however,  we  will  only  show  two  examples  of 
reconstruction  where  the  number  of  crossings  on  the  sampling  lines  are  smaller  than  the  number 
of  Fourier  coefficients  of  the  one-dimensional  signal  associated  with  the  sampling  line.  Figure 
(3.17a)  shows  the  reconstructed  version  of  the  eye  picture  from  4,6  and  8  thresholds,  via 
the  iterative  algorithm  for  M  =  64, 128.  Similar  to  the  previous  example,  we  chose  32 

equidistant  sampling  lines,  and  20  iterations  for  each  one-dimensional  reconstruction.  The 
relaxation  parameters  were  chosen  to  be  \\  =  1,  A2  =  1.75.  Our  second  example  shown 
in  figure  (3  17b)  deals  with  the  reconstruction  of  the  63  x  63  cman  picture  from  4,6  and  8 
thresholds  with  M  =  128,256.  We  used  64  equidistant  sampling  lines  of  unit  slope,  the 

number  of  iterations  for  each  line  was  20,  and  the  relaxation  parameters  were  chosen  to  be 
Ai  =  1,  A2  =  1.75. 

As  seen  in  both  of  the  above  examples,  increasing  the  number  of  thresholds  affects  the  quality 
of  reconstruction  in  a  more  substantial  way  than  increasing  M.  We  have  found  experimentally 
that  if  the  sixth  step  of  the  algorithm  dealing  with  projection  onto  the  space  domain  is  modified 


Figure  3.17:  Reconstructed  version  of  the  eye  and  cman  pictures  via  the  iterative  algorithm. 
The  numbers  of  thresholds  are  4,6  and  8,  increasing  from  left  to  right.  The  number  of  equally 
spaced  samples  on  the  lines  increases  from  top  to  bottom,  and  are  64,  128  for  the  eye  picture, 
and  128  and  256  for  the  cman  picture. 
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to  accelerate  convergence.  Recall  that  the  sixth  step  of  the  algorithm  projects  the  sequence 
g  onto  C2,  the  set  of  M  point  real  sequences  satisfying  the  space  domain  constraint.  This 
implies  that  if  the  value  of  the  unknown  function  is  constrained  to  be  smaller  than  tu(n)  and 
larger  than  tj(n),  and  is  outside  the  interval  [t|(n),  tu(n)],  its  projected  value  is  either 

tu(n)  or  ti(n),  depending  on  the  value  of  g(xf)-  Although  the  the  scheme  shown  in  equation 
(3.7)  is  not  a  projection  and  although  the  mean  square  error  does  not  always  decrease  from  one 
iteration  to  the  next,  by  “projecting”  onto  instead  of  tt(n)  or  tu(n),  we  can  somewhat 

speed  up  the  convergence.  This  can  be  explained  by  formulating  the  problem  in  a  statistical 
framework  [49],  but  we  shall  not  go  over  it  here.  The  important  point  here  is  that,  although  in 
most  reconstructions,  carried  out  via  the  above  modified  version  the  algorithm,  converged  to  a 
sequence  satisfying  both  the  space  and  frequency  domain  constraints,  the  mean  square  error  of 
the  solution  was  almost  identical  to  the  solution  obtained  via  the  over  relaxation  method.  This 
suggests  that  the  reason  behind  the  poor  quality  of  the  reconstruction  for  a  small  number  of 
thresholds  is  not  the  fact  that  the  algorithm  did  not  converge,  but  that  there  are  many  signals 
satisfying  the  space  and  frequency  domain  constraints;  the  algorithm  has  simply  converged  to 
one  of  them.  Increasing  the  number  of  thresholds  or  increasing  M,  imposes  more  constraints 
on  the  signal,  thus  reducing  the  number  of  sequences  satisfying  both  the  space  and  frequency 
domain  requirements.  As  we  said  earlier,  if  the  number  of  crossings  on  a  sampling  line  exceeds 
the  number  of  Fourier  harmonics  of  the  one-dimensional  signal  associated  with  it,  the  solution, 
obtained  via  the  iterative  algorithm,  becomes  unique  as  M  — *  00. 

Before  we  end  this  section,  let  us  briefly  compare  the  iterative  algorithm  with  other  recon¬ 
struction  algorithms  proposed  in  this  chapter.  Unlike  the  iterative  approach  which  needs  all 
the  intersections  of  sampling  lines  with  the  crossings  of  the  signal  in  order  to  derive  the  space 


domain  constraint  correctly,  all  of  our  other  algorithms  can  use  a  subset  of  the  intersections 
for  successful  reconstruction.  Clearly,  the  quality  of  reconstruction  via  the  iterative  algorithm 
is  not  as  high  as  the  least-squares  approach  or  the  non-iterative  line-by-line  reconstruction; 
However,  one  must  bear  in  mind  that  reconstructions  using  the  latter  techniques  use,  at  least 
in  principle,  the  exact  location  of  the  crossings,  whereas  the  iterative  algorithm,  in  effect,  uses 
the  rather  coarsely  quantized  version  of  the  crossings  as  long  as  M  is  finite.  We  will  study  these 
quantization  issues  in  more  depth  in  Chapter  5.  Finally,  the  iterative  algorithm  is  free  of  many 
limitations  of  the  other  schemes  proposed  earlier  in  this  chapter.  More  specifically, 

•  Its  storage  requirement  is  not  as  stringent  as  the  QR  decomposition,  thus,  enabling  us  to 
process  much  larger  images.  In  fact,  it  is  the  only  algorithm  in  this  chapter  which  can 
reconstruct  large  images  from  multiple  level  crossings  (as  opposed  to  sinusoid  crossings). 

•  It  is  not  as  computationally  intensive  as  the  QR  decomposition,  conjugate  gradient  or  the 
non-iterative  line-by-line  reconstruction  scheme. 

•  It  is  considerably  more  robust  than  the  recursive  algorithm  of  section  (3.2),  which  is  only 
capable  of  reconstruction  from  sinusoid  crossings. 

•  It  is  more  robust  than  the  other  algorithms  in  the  sense  that  if  the  number  of  thresholds  is 
not  large  enough  to  result  in  a  sufficient  number  of  crossings,  it  can  still  carry  out  approx¬ 
imate  reconstruction.  Indeed,  this  graceful  degradation  in  the  quality  of  reconstruction  is 
a  feat  jre  which  is  lacking  in  all  the  other  algorithms  we  have  discussed  so  far. 

In  Chapter  5,  we  will  discuss  the  quantization  characteristics  of  the  iterative  algorithm  and  the 
linear  least-squares  approach  in  detail. 


3.4  Discussion 


In  this  chapter,  we  proposed  three  major  reconstruction  algorithms  for  the  semi-implicit 
sampling  strategy  of  Chapter  2.  Our  first  approach  involved  solving  a  linear  system  of  equa¬ 
tions  via  QR  decomposition  or  the  conjugate  gradient  algorithm.  Reconstruction  via  QR  de¬ 
composition  is  more  stable  than  all  of  our  other  algorithms.  However,  its  storage  requirements 
are  large  enough  to  limit  its  applicability  to  reconstruction  problems  with  a  large  number  of 
Fourier  coefficients.  In  fact,  for  a  signal  with  N  x  N  region  of  support  in  the  Fourier  domain, 
the  storage  requirement  of  the  QR  decomposition  is  proportional  to  N6,  and  its  computation 
requirements  are  proportional  to  N6.  The  conjugate  gradient  algorithm  does  not  need  to  store 
the  entire  linear  least-squares  matrix.  It  is,  therefore,  less  storage  intensive.  However,  its  major 
drawback  is  that  unless  the  reconstruction  problem  is  extremely  well  conditioned,  it  does  not 
result  in  satisfactory  reconstruction.  Thus,  we  have  only  been  able  to  use  it  for  recovery  of 
signals  from  their  sinusoid  crossings. 

The  recursive  algorithm  proposed  in  section  (3.2)  takes  advantage  of  the  geometry  of  the 
interpolation  points  obtained  via  the  semi-implicit  sampling  strategy.  It  is  not  as  general  as 
the  linear  least-squares  approach,  since  it  only  deals  with  sampling  situations  in  which  all 
the  sampling  lines  have  fixed  integer  slopes.  Furthermore,  since  it  computes  the  coefficients 
recursively,  small  errors  in  the  coefficients  obtained  in  the  initial  stages  of  the  algorithm  are 
propagated  and  magnified  in  the  final  steps,  resulting  in  instability.  Although  it  is  somewhat 
more  stable  than  the  conjugate  gradient  algorithm,  it  can  only  be  applied  to  well  conditioned 
problems  such  as  reconstruction  from  sinusoid  crossings. 

The  non-iterative  line-by-line  reconstruction  scheme  of  section  (3.3)  is  more  stable  than  the 
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conjugate  gradient  or  the  recursive  algorithm,  and  its  computation  and  storage  requirements 
are  not  as  demanding  as  the  QR  decomposition.  Its  major  drawback  however,  is  that  the 
number  of  crossings  it  needs  for  reconstruction  exceeds  the  theoretical  requirements  of  Chapter 
2.  Thus  for  reconstruction  from  multiple  level  crossings,  the  number  of  thresholds  required  by 
this  algorithm  is  generally  larger  than  our  other  algorithms. 

Finally,  the  iterative  algorithm  of  section  (3.3.1)  shares  many  of  the  characteristics  of  the 
non-iterative  line-by-line  reconstruction  algorithm  except  that  it  is  capable  of  approximate 
reconstruction  in  situations  where  the  number  of  crossings  do  not  satisfy  the  theoretical  re¬ 
quirements  of  Chapter  (2).  However,  it  requires  all  the  intersections  of  sampling  lines  with  the 
crossings  of  the  signal  in  order  to  derive  the  space  domain  constraint  correctly. 

In  spite  of  the  variety  of  reconstruction  algorithms  proposed  in  this  chapter,  the  major 
drawback  of  the  semi-implicit  sampling  scheme  is  the  fact  that,  in  general,  we  are  never  guar¬ 
anteed  to  get  enough  intersections  between  the  sampling  lines  and  the  level  crossings  in  order 
to  satisfy  the  minimum  theoretical  requirements  of  the  Chapter  2.  Of  course,  we  can  propose 
various  guidelines  in  order  to  increase  the  number  of  intersections  as  much  as  possible.  For  in¬ 
stance,  one  way  would  be  to  choose  the  slopes  of  our  sampling  lines  perpendicular  to  the  general 
direction  of  the  crossing  contours.  More  specifically,  if  the  harmonic  content  of  a  :eal  signal  at 
spatial  frequency  (p,  q)  is  very  strong,  we  would  expect  the  signal  to  vary  with  frequency  (p,  q) 
along  lines  of  slope  *  and  thus  choosing  the  sampling  lines  with  slope  *  would  result  in  large 
number  of  intersections. 

Another  way  to  increase  the  likelihood  of  satisfying  the  theoretical  requirements  is  to  choose 
sampling  lines  which  result  in  the  least  stringent  sampling  requirements.  For  instance,  in 
situations  where  all  the  sampling  lines  have  equal  slope  j£,  the  ratio  between  the  maximum 
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number  of  points  on  a  sampling  line  and  its  length  is  given  by: 


2  _  imm  _+_n)  +_1J2  2  _ _ ' 

Pmaz  n2  +  m2  L  1  +  (n/m)2  J 


Thus  for  jjj  >  1,  larger  slopes  result  in  less  stringent  sampling  requirements  whereas  for  0  < 
~  <  1,  smaller  slopes  result  in  less  stringent  sampling  requirements. 

Clearly,  both  of  the  above  guidelines  are  somewhat  qualitative  and  by  no  means  do  they 
guarantee  that  the  number  of  intersections  of  the  sampling  lines  with  the  level  crossings  sat¬ 
isfies  the  theoretical  requirements  of  Chapter  2.  The  major  drawback  of  the  semi-implicit 
approach  still  remains,  since  it  is  incapable  of  providing  us  with  sufficient  conditions  for  unique 
reconstruction  from  an  arbitrarily  small  number  of  thresholds.  Indeed,  this  has  been  our  main 
motivation  for  deriving  the  implicit  sampling  scheme  of  the  next  chapter. 
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Chapter  4 

Implicit  Approach  to 
Reconstruction 

In  Chapter  2,  we  proposed  a  variety  of  semi  implicit  sampling  strategies  for  the  reconstruc¬ 
tion  of  multidimensional  signals  from  their  crossings  with  arbitrary  functions.  More  specifically, 
the  intersection  of  sampling  lines  of  rational  slope  with  crossing  contours  was  used  to  interpo¬ 
late  the  bivariate  polynomial  associated  with  the  signal.  As  we  mentioned  earlier,  the  major 
drawback  of  the  line  sampling  strategy  for  reconstruction  from  function  crossings  is  the  fact 
that,  in  general,  one  is  never  guaranteed  to  get  enough  intersections  between  the  sampling  lines 
and  signal  crossings  to  satisfy  the  conditions  of  Corollaries  (2.1),  (2.2)  and  (2.3).  As  we  saw 
in  section  (3.4),  we  can  improve  the  likelihood  of  satisfying  the  sampling  requirements  of  these 
theoretical  results  by  careful  choice  of  the  slope  and  position  of  sampling  lines.  In  addition,  in 
section  (3.4),  we  found  that  under  certain  conditions,  the  number  of  intersections  of  sampling 
lines  with  certain  functions  such  as  sinusoids,  becomes  predictable.  In  spite  of  these  observa¬ 
tions,  the  basic  problem  still  exists;  that  is,  for  arbitrary  functions,  the  intersections  of  sampling 
lines  and  level  crossings  might  not  result  in  a  sufficient  number  of  interpolation  points  to  satisfy 
our  theoretical  results  of  Chapter  2.  Therefore,  it  seems  natural  to  search  for  more  general 
strategies  for  reconstruction  of  multidimensional  signals  from  their  crossings  with  an  arbitrary 


number  of  thresholds.  The  derivation  of  such  a  scheme  will  be  our  major  goal  in  this  chapter. 

As  we  mentioned  in  Chapter  2,  the  inherent  difficulty  in  bivariate  interpolation  is  the  lack 
of  existence  of  Chebychef  systems  in  R 2.  Our  approach  in  Chapter  2  was  to  constrain  the 
location  of  our  interpolation  points.  Another  approach,  which  will  be  the  basis  for  our  implicit 
sampling  strategy  of  this  chapter,  consists  of  conditionally  regular  interpolation  methods.  An 
interpolation  method  is  called  regular  if  it  is  uniquely  solvable  for  any  selection  of  interpolation 
points.  Conditionally  regular  interpolation  methods  are  solvable  not  for  all  selection  of  points, 
but  only  for  most  of  them  [31].  Roughly  speaking,  they  are  uniquely  solvable  with  probability 
one.  For  methods  of  this  type,  if  one  has  a  concrete  problem,  and  selects  the  interpolation 
points  at  random,  it  will  be  extremely  likely  that  the  problem  will  be  solvable. 

In  section  (4.1),  we  will  take  the  above  approach  to  derive  sufficient  conditions  under  which 
the  problem  of  reconstruction  from  level  crossings  is  almost  always  uniquely  solvable.  Section 
(4.2),  describes  two  reconstruction  algorithms  for  the  results  of  section  (4.1).  It  is  worthwhile 
to  mention  that  while  the  results  of  Chapter  2  are  applicable  to  the  more  general  problem  of 
reconstruction  from  crossings  with  arbitrary  functions,  the  results  of  this  section  are  specifically 
tailored  to  reconstruction  from  level  crossings  or  crossings  with  any  pre-specified  function  that 
lies  in  the  band  of  the  signal  under  consideration.  On  the  other  hand,  unlike  the  semi-implicit 
approach,  the  implicit  approach  of  this  chapter  is  capable  of  reconstruction  from  an  arbitrary 
number  of  level  crossings. 

4.1  Theoretical  Results 

In  this  section,  we  apply  algebraic  geometric  concepts  to  derive  conditions  under  which 
multidimensional  signals  can  be  reconstructed  from  their  level  crossings  with  an  arbitrary  num- 
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ber  of  thresholds.  More  specifically,  we  will  show  that  under  mild  conditions,  a  signal  with 
(2 N  +  1)  x  (21V  +  1)  region  of  support  in  the  Fourier  domain  can  be  uniquely  reconstructed 
from  almost  any  k  points  from  its  a  level  crossings,  and  (21V  +  l)2  -  k  points  from  its  0  level 
crossings.  As  we  will  see,  the  extension  of  the  above  result  to  the  problem  of  reconstruction  from 
more  than  two  threshold  crossings  is  rather  straightforward.  Furthermore,  we  can  extend  our 
basic  result  to  the  problem  of  reconstruction  from  sinusoid  crossings,  provided  the  frequency  of 
the  sinusoid  under  consideration  lies  within  the  bandwidth  of  our  signal. 

Consider  a  two-dimensional,  real,  BLP  signal,  /( i,  y)  with  periods  Tx  and  Tv  in  x  and  y 
directions  respectively.  The  Fourier  series  representation  of  /(x,  y)  is  given  by: 

f(x,y)  =  £  £  F(kt,k2)  ej2r{#  +  (4.1) 

For  simplicity  we  will  assume  F{k\,  k2)  to  have  a  square  region  of  support  of  size  (21V  +  1)  x 
(21V  +  1),  and  Tz  and  T„,  the  periods  in  x  and  y  directions  to  be  equal  to  1.  Our  results  can 
be  easily  modified  for  cases  when  Tx  and  Tv  are  not  equal  to  1  and  F{k\,  k2)  has  a  rectangular 
region  of  support. 

Since  /(x,  y)  is  real,  its  Fourier  series  coefficients  are  conjugate  symmetric.  That  is: 


F{kuk2)  =  F*(-fci, -fcj) 


l*il»l*il  <  N 


(4.2) 


Using  this  relationship,  we  can  rewrite  our  signal  as  a  trigonometric  function  with  real  coeffi¬ 
cients  in  the  following  manner: 


/(*,*)  =  ^(0,0)]  + 


(4.3) 


N 


2  <  i?e[F(0,A:2)]cos(2jrfc2y)  -  /m[F(0,  fc2)]st'n(2jrfc2y)  + 

U»=i 

N  N 

X]  H  rte[F(fc],fc2)]coa|2*(fciz  +  k2y)}  -  Im\F{ki,k2)}8xn{2^(kix  +  fc2y)] 
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as  opposed  to  square. 

As  we  mentioned  earlier,  our  main  goal  is  to  show  that  almost  any  0  <  k  <  (2 N  +  l)2 
points  from  the  set  of  a  level  crossings  and  (2N  +  l)2  -  k  points  from  the  0  level  crossings 
are  sufficient  for  unique  specification  of  f(x,  y)  provided  a  ^  0.  In  order  to  express  this  in  a 
mathematical  framework,  let  us  define  the  sets  Aa  and  A$  as  a  and  0  level  crossings  of  /,  in 
the  following  manner: 

^a(C)  =  {(x, y)  €  C2  |  /(x,y)  =  a} 

Aa(R)  =  {(z,y)  €  R2  |  /(*,»)  =  a} 

=  Aa(C)f)R2  (4.5) 

Af)(c)  =  {{x,y)eC2  |  /(*, y)  =  0 } 

Ap{R)  =  {(x,y)€  R2  |  /(*,y)  =  d} 

=  Afi{C)  n  R2  (4.6) 

Furthermore,  let  B{R)  C  R2(2N+1)2  denote  the  cross  product 

B(R)  =  Aka(R)  x  A[02N+l)2-k(R) 

where  A*(R)  is  k  times  the  cross  product  of  Aa(f?)  with  itself,  and  A (2N+l'1  k(R)  is  (2JV+l)2-fc 
times  the  cross  product  of  Ap(R)  with  itself. 

Having  made  the  above  definitions,  we  can  now  state  our  goal  as  finding  the  conditions 
under  which  real  zeros  of  the  trigonometric  polynomial  p(xi,  yi,  ...,X(2jv+i)3>  S/(2Ar+ip)  become 
of  measure  zero  in  the  set  B(R)  C  ft2(2^+i)2  To  formulate  our  problem  in  an  algebraic 
geometric  framework,  we  must  change  some  variables  to  express  our  trigonometric  polynomials 
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as  algebraic  ones.  This  is  primarily  because  algebraic  geometry  only  deals  with  zeros  of  algebraic 
polynomials  rather  than  trigonometric  ones.  To  this  end,  let 

stn{2jrx)  =  w 

cos(  2jtx)  =  2 
sin(2iry)  =  ti 
cos(2?ry)  =  t; 

so  that  the  trigonometric  polynomials  /(x,  y)  and  p(xi,t/i, ...» X(2/v+i)2,  y{2N  +  i)*)  can  be  rewrit¬ 
ten  as: 

/(*.»)  =  /(u>,r,u,t>) 

p(*l»  Vl>  •>a:(2tf+l)a>y(2JV-l-l)s)  =  p(u'lizl>ul>t,l>  *(2N  +  1)2>  w(2JV  +  1)2>  v(2JV  +  I)2) 

With  these  new  variables,  the  sets  Aa  and  Ap  can  be  rewritten  in  the  following  way. 


*«(C)  = 


(w,z,u,v)  e  C4  I  /(w,z,u,v)  =  a 
w2  +  z2  —  1 
u2  -I-  V2  =  1 


Aa{R)  =  Aa(C)  HR4 
(w,z,u,v)  €  C4  I  f(vu,z,u,v)  =  0 
MC)  =  {  u}2  4-  z2  —  1 

u2  -I-  v2  =  1 

Aff{R)  =  ap{c)  n  R* 

B(C)  =  A*(C)  X  A<2Ar+,)’-*(C) 
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(4.7) 
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B{R)  =  B{C)  nR2 t2W+1)s 

=  Ak(R)  x  A[™+l]2-k(R) 

Since  there  is  a  one  to  one  correspondence  between  elements  of  the  sets  Aa  and  Aa,  Ap 
and  Ap,  and  B  and  B,  our  major  goal,  which  was  to  find  conditions  under  which  real  zeros 
of  p(ii,yx,  •■•,£(2jv+i):Ji  y(2Jv+i)2)  have  measure  zero  in  B(R ),  becomes  equivalent  to  finding 
conditions  under  which  real  zeros  of  p(w i,  zj,  •••,  u^n+i)2’  v(2Ar+i)2)  have  measure  zero  in  B(R). 
The  major  implication  of  changing  variables  is  the  fact  that,  unlike  Aa,  Ap  and  B ,  the  sets 
Aa,  Ap  and  B  are  algebraic  sets.  Algebraic  sets  are  important  elements  in  algebraic  geometry, 


and  we  will  take  full  advantage  of  their  properties  in  proving  the  major  result  of  this  chapter. 
Thus,  before  we  state  our  major  theorem,  it  is  worthwhile  to  go  over  few  definitions  in  algebraic 
geometry  [50]. 

Definition  4.1  Let  1C  be  any  field  and  An(K)  or  simply  An  denote  the  Cartesian  product  of  1C 
with  itself  n  times.  That  is  An(K)  is  the  set  of  n-  tuples  of  elements  of  K.  If  the  coefficients 
of  a  polynomial  F  are  in  1C ,  i.e.  F  E  IC\X\,  ...,Xn],  a  point  P  =  (aj ,  ...,an)  £  An(IC)  is  a  zero 
of  F  if  F[P)  =  F(ai,  ...,an)  =  0.  If  S  is  any  set  of  polynomials  in  K[Xi,...,Xn\,  we  define 
V (S)  to  be  the  set  of  common  zeros  of  all  the  polynomials  in  S.  That  is: 


V{S)  =  {PeAn\F{P)  =  0,  VFeS} 

A  subset  X  C  An(IC)  is  an  affine  algebraic  set  or  simply  an  algebraic  set,  if  X  —  V(S)  for 
some  S. 

Thus  Aa(R)  and  Ap(R)  are  algebraic  sets  in  R*,  and  Aa{C )  and  A/j(C)  are  algebraic  sets 
in  C* .  An  algebraic  set  may  be  the  union  of  several  smaller  algebraic  sets.  An  algebraic  set 
V  C  An  is  reducible  if 

V  =  Vx  u  V2 

where  Vj  and  Vi  are  algebraic  sets  in  A"  and 

Vi*V  i  =  1,2 
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V\  and  V2  are  referred  to  as  components  of  V.  If  V  is  not  reducible,  it  is  irreducible.  It  can  be 


1. 

2. 


shown  that  an  algebraic  set  is  the  union  of  a  finite  number  of  irreducible  algebraic  sets.  More 
specifically,  we  have: 


Theorem  4.1  f50]:  Let  V  be  an  algebraic  set  in  An(K).  Then  there  are  unique  irreducible 
algebraic  sets  V1,...,Vm  such  that 

V  =  V,  U  ...  U 


and 

Vi  <t  Vj,  Vi 


The  proof  is  straightforward  and  is  included  in  [50j .  An  important  consequence  of  the 
above  theorem  is  the  fact  that  given  an  irreducible  algebraic  set  V  in  An(/C),  the  zeros  of  the 
polynomial  F  e  K\Xi,  ...,X„]  are  either  of  measure  zero  in  V,  or  contain  all  the  points  in  V. 
As  we  will  see,  this  fact  will  be  crucial  in  proving  our  major  theorem  which  can  be  stated  in 
the  following  manner: 


Theorem  4.2  The  set  of  real  zeros  of  the  polynomial  p(w\, ...,  v^jv+ip)  *8  °f  measure  zero  in 
the  set  B{R)  provided  the  following  conditions  hold: 

Aa(R)  and  Ap{R)  have  maximal  topological  dimensions. 

The  polynomials 

2N  2  N 

ga{WuW2)  =  £  £  F^h- N,k2- N)W['Wp 

*1  =  0  *2=0 


2  N  2  N 

9fi(WltW2)  =  J2  £  *>(*1 
*,=0  *2=0 


with 


.*»)  =  | 
.*>)  =  { 


F{ 0,0)  -  a 

F{kr,ki) 

F{ 0,0) -0 

F(kt,k2) 


are  irreducible  over  the  set  of  complex  numbers. 


N,  k2  -  N)W*lW$* 

k\  =  k2  =  0 
elsewhere 

ki  —  k2  —  0 
elsewhere 
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The  outline  of  the  proof  is  as  follows:  We  will  first  use  the  first  condition  to  show  that  B[R)  is  an 
irreducible  algebraic  set  over  reals.  Then  either  zeros  of  p  are  of  measure  zero  in  B(R)  or  all  the  points 
in  B(R)  are  zeros  of  p.  Since  our  objective  is  to  prove  that  real  zeros  of  p  are  of  measure  zero  in  B(R), 
all  we  have  to  show  is  that  there  is  at  least  one  point  in  B[R)  at  which  p  does  not  vanish.  As  we  will 
see,  the  second  condition  of  the  theorem  will  be  used  to  show  this. 

To  begin,  notice  that  the  polynomials  ga(W1,W2)  and  gp(Wlt  W2)  are  related  to  our  BLP  signal 
}(x ,  y)  via  the  change  of  variables 

Wi  =  e]2nx 
W2  =  ej2,rw 

in  the  following  manner: 

/(*,»)  ~  «  =  W~N  W2N  w2) 
fix,  y)-0  =  W~N  W~N  gp(WuW2) 

Thus,  there  is  a  one  to  one  correspondence  between  zeros  of  yu(W i,W2)  (resp.  gp(Wi,W2))  and  Aa(C) 
(resp.  Ap(C).  Therefore,  the  second  condition  of  the  theorem  implies  that  A„  and  Ap  are  also  irre¬ 
ducible.  Considering  equation  (4.7),  since  B(C)  is  the  Cartesian  product  of  A„(C)s  and  A^(C)s,  we  can 
conclude  that  B(C )  is  also  irreducible  over  complex  numbers. 

Similarly,  since  Au(fi)  and  Ap(R)  have  maximal  topological  dimensions,  so  does  B(R).  Therefore 
considering  Theorem  (C.l)  of  Appendix  (C),  and  taking  into  account  that  B(C)  is  irreducible,  we  can 
conclude  that  B{R)  is  Zariski  dense  in  B(C).  This  means  that  every  polynomial  that  vanishes  on  B(R) 
vanishes  on  ail  of  B(C).  This,  together  with  the  fact  that  B(C)  is  irreducible,  implies  that  B(R)  is 
irreducible  over  complex  numbers  and  thus  the  reals.  Therefore,  in  order  to  show  that  the  real  zeros  of 
p(u>i,  ...,i)(w+l|j)  are  of  measure  zero  in  B(R),  we  merely  have  to  show  that  there  exists  at  least  one 
point  in  B(R)  at  which  p  does  not  vanish.  That  is,  p  does  not  vanish  identically  on  B(R). 

Since,  by  hypothesis,  ya(W j,  W2)  and  gp{W^,  W2)  are  irreducible  over  complex  numbers,  taking  into 
account  modified  version  of  Bezout’s  theorem  of  Appendix  (B),  we  can  conclude  that  any  81V2  + 1  samples 
of  Aa  will  enable  us  to  specify  f(x,  y)— a  to  within  a  scale  factor  [51].  This  means  the  (81V2 -t-l)  x  (2JV+1)2 
matrix  associated  with  any  81V2  +  1  points  of  A„,  Da  has  rank  (21V  +  l)2  —  1,  and  its  null  vector  is 
specified  by  the  coefficients  of  /( i,  y)  —  a.  Similarly,  Dp,  the  81V2  +  1  x  (21V  +  l)2  matrix  associated 
with  any  81V2  +  1  points  of  Ap  has  rank  (2 IV  +  l)2  —  1,  and  its  null  vector  is  given  by  the  coefficients 
of  /( x,  y)  —  0.  Since  by  hypothesis  a  ^  0,  the  direction  of  the  null  vectors  of  D„  and  Dp  are  different 
from  each  other.  Therefore,  there  exists  at  least  one  combination  of  k  rows  from  D„  and  (21V  +  l)2  —  k 
rows  from  Dp  which  result  in  a  full  rank  (21V  -I-  l)2  x  (21V  -t-  l)2  matrix  with  non  zero  determinant  p. 
Hence  p  does  not  vanish  identically  on  B(R). 

Since  B(R)  is  irreducible  over  reals  and  p  does  not  vanish  on  it  identically,  the  real  zeros  of  p  must 
be  of  measure  zero  in  B(R).  QED. 

□ 

The  above  theorem  implies  that  two-dimensional  signals  with  (21V  +  1)  x  (21V  +  1)  region 
of  support  in  the  Fourier  domain  can  be  uniquely  specified  from  almost  any  0  <  k  <  (2 N  +  l)2 
points  from  their  a  crossings  and  (2rV  -r  l)2  -  k  points  from  their  crossings,  provided  they 
satisfy  the  conditions  of  the  theorem.  Thus,  it  is  natural  to  question  the  strictness  of  these 


conditions  in  practical  situations.  The  first  condition  of  the  theorem  requires  Aa{R)  and  Ap(R) 
to  have  maximal  topological  dimensions.  Since  the  complex  topological  dimension  of  Aa(C) 
(resp.  ^(C))  is  1,  and  its  real  topological  dimension  is  2,  maximal  dimension  for  Aa(R)  (resp. 
Ap(R))  is  1.  This  implies  that  the  points  (*,  y)  satisfying  f(x,y)  =  a  (  resp.  /(x, y)  =  0) 
must  contain  at  least  a  curve  and  can  not  simply  be  isolated  points  in  the  x  -  y  plane.  In 
practice,  this  condition  can  be  easily  satisfied  by  choosing  our  threshold  a  (resp.  /?)  in  such  a 
way  that  the  corresponding  threshold  contours  form  curves  and  not  isolated  points. 

The  second  condition  of  Theorem  (4.2)  is  also  easily  satisfied  in  practice.  Since  throughout 
our  derivation  we  assumed  /(x,y)  to  be  real  and  F’(fci,Jbj)  to  be  conjugate  symmetric,  our 
approach  is  to  show  that  the  set  of  reducible  polynomials  is  of  measure  zero  in  the  space  of  all 
the  bivariate  polynomials  of  the  form: 

2  N  2  N 

p{w,z)  -  52  »(«./) 

t=0  j — 0 

a"(«,  j)  =  a{2N-i,2N-j) 

As  it  turns  out,  this  is  a  simple  extension  of  the  fact  that  the  set  of  reducible  polynomials  are 
of  measure  zero  in  the  space  of  all  the  bivariate  polynomials  with  complex  coefficients  [52,53]. 

Theorem  (4.2)  can  be  easily  extended  to  recovery  of  multidimensional  signals  from  more 
than  two  threshold  crossings.  More  specifically,  we  can  show  that  given  m  distinct  thresholds, 
ti,  almost  any  distribution  of  (2 N  + 1)2  points  among  the  thresholds  will  result  in  unique 

reconstruction  of  the  signal  under  consideration  provided 

1.  Not  all  the  points  are  chosen  on  one  threshold. 

2.  For  1  <  i  <  m,  the  set  of  points  satisfying 


has  maximal  topological  dimension  (i.e.  consists  of  curves  as  opposed  to  isolated  points). 

3.  For  1  <  i  <  m,  the  polynomial  associated  with 

/(*,  y)  =  U 

is  irreducible  over  complex  numbers. 

Furthermore,  Theorem  (4.2)  implies  that  almost  any  {2N  + 1)2  samples  of  one  level  crossings 
of  a  signal  with  (2 N  +  1)  x  (2N  +  1)  region  of  support  in  the  Fourier  domain  is  sufficient  for 
its  unique  reconstruction  to  within  a  scale  factor,  provided  the  polynomial  associated  with  the 
signal  is  irreducible.  This  is  in  contrast  with  (16JV2  +  1)  samples  required  by  the  theoretical 
results  of  Curtis  [51],  or  8 iV2  +  1  samples  required  by  the  application  of  the  modified  form  of 
Bezout’s  theorem  as  in  [39], 

Finally,  it  is  not  hard  to  see  that  the  basic  idea  in  Theorem  (4.2)  can  be  extended  to 
reconstruction  from  crossings  with  BLP  functions,  whose  bandwidth  lie  within  the  bandwidth 
of  the  signal.  Examples  of  reconstruction  via  the  implicit  sampling  scheme,  described  in  this 
section,  are  included  in  the  next  section. 

4.2  Reconstruction  Algorithms 

In  the  previous  section,  we  found  that  most  multidimensional  signals  can  be  uniquely  re¬ 
constructed  from  samples  on  contours  corresponding  to  two  or  more  level  crossings.  Unlike  the 
semi-implicit  line  sampling  strategy  of  Chapter  2,  the  sampling  scheme  described  in  section 
(4.1)  is  truly  implicit.  In  other  words,  whereas  the  interpolation  points  of  the  semi-implicit 
approach  must  lie  on  lines  of  rational  slope,  for  the  implicit  approach,  there  is  no  pre-imposed 
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pattern  or  structure  to  the  samples. 


V 


In  this  section,  we  will  propose  two  methods  to  reconstruct  multidimensional  signals  from 
their  level  crossings,  via  the  implicit  sampling  approach  of  the  previous  section.  The  most 
straightforward  way,  is  to  solve  a  linear  (possibly  overdetermined)  system  of  equations  of  the 
form: 

N  N 

/(*,,  y.)  =  £  £  F(k,,k2y2*lk>*’  +  *»*•)  =  t,  1  <  *  <  M,  1<J<P 

fc,  =  -N 

where  (i,,  y,)  correspond  to  the  threshold  values  t;.  Theoretical  results  of  the  previous  section 
guarantee  the  uniqueness  of  the  solution  provided  the  signal  under  consideration  satisfies  the 
mild  conditions  of  Theorem  (4.2).  An  example  of  reconstruction  of  the  original,  eye  picture 
of  figure  (2.10)  from  samples  of  its  two-level  crossings  is  shown  in  figures  (4.1)  and  (4.2). 
Figure  (4.1)  shows  the  two-level,  non-uniform,  amplitude  quantized  version  of  the  eye  picture 
and  the  crossing  contours  associated  with  thresholds  75  and  145.  Since  the  eye  picture  has 
961  coefficients,  invoking  Theorem  (4.2),  almost  any  k  points  from  its  75  level  crossings  and 
961  -  k  samples  from  its  145  level  crossings  are  sufficient  for  its  unique  reconstruction.  The 
reconstructed  version  of  the  eye  picture  using  the  linear  least-squares  approach  is  shown  in  figure 
(4.2).  QR  decomposition  was  used  with  526  samples  on  the  threshold  contour  corresponding  to 
75,  and  426  samples  on  the  145  contour.  Although  other  linear  least-squares  algorithms  such  as 
conjugate  gradient  algorithm  need  less  storage  than  QR  decomposition,  they  are  generally  not 
stable  enough  for  the  problem  of  reconstruction  from  multiple  level  threshold  crossings,  unless 
the  number  of  thresholds  is  very  large.  Our  second  approach  to  reconstruction  is  similar  to 
the  iterative  algorithm  of  section  (3.3.1).  The  remaining  part  of  this  section  is  devoted  to  this 
algorithm. 


U! 


104 


Figure  4.1:  (a)  2  level  non  uniform  amplitude  quantized  version  of  the  original  eye  picture;  (b) 
Threshold  contours  associated  with  the  2  level  crossings. 


Figure  4.2:  Reconstructed  version  of  the  eye  picture  from  its  level  crossings  at  75  and  145  via 
the  QR  decomposition. 


4.2.1  The  Iterative  Approach 


Our  second  approach  to  the  problem  of  reconstruction  from  implicit  samples  of  level  cross¬ 
ings  is  iterative  and  is  based  on  Theorem  (3.1)  of  section  (3.3.1).  Intuitively,  this  strategy  is 
based  on  imposing  the  space  and  frequency  domain  constraints  successively.  The  space  domain 
constraint  is  derived  from  the  quantized  threshold  contours;  the  frequency  domain  constraint 
is  due  to  the  bandlimitedness  of  the  signal  under  consideration.  We  will  now  describe  the  algo¬ 
rithm  for  reconstructing  a  BLP  signal  with  (2 N  +  1)  x  (2 N  +1)  region  of  support  in  the  Fourier 
domain.  We  assume  that  all  the  crossing  contours  of  the  signal  associated  with  p  thresholds 
*1  <  *2  <  •••  <  *p  are  quantized  on  an  M  x  M  grid  where  M  >  2 N  +  1,  and  that  the 
intensity  of  the  signal  lies  in  the  range  [fo,*p+i]-  The  steps  of  the  algorithm  are  as  follows: 

1.  Deduce  the  range  of  intensity  of  the  signal  on  the  nodes  of  the  Af  x  M  grid  using  the 
position  quantized  threshold  crossings.  The  quantization  process  is  shown  pictorially  in 
figure  (4.3).  As  shown,  quantized  threshold  contours  trace  outside  of  the  boundary  of 
level  crossings  on  the  nodes  of  the  grid.  Thus,  as  long  as  there  are  a  minimum  of  two 
contours  corresponding  to  two  different  thresholds,  the  intensity  range  for  the  nodes  on 
the  MxM  grid  can  be  found1.  Therefore,  if  ti(ni,n2)  and  tu(n\,n2)  denote  the  lower 
and  upper  bound  for  the  (nj,n2)th  point  on  the  grid,  we  have: 


I 


ti{ni,n2)  <  ^  tu{ni,n2),  n1,n2  =  0,  ...,  M-l 


It  is  worthwhile  to  mention  that  the  quantized  threshold  contours  can  also  be  chosen  to 
be  inside  of  the  actual  contours.  In  general,  our  quantization  scheme  must  be  chosen  so 


‘The  fact  that  we  need  at  least  two  thresholds  Is  a  direct  consequence  of  Theorem  (4.2).  Intuitively,  if  we  only 
have  contours  associated  with  one  threshold,  there  is  no  way  to  determine  which  parts  of  the  signal  are  smaller 
than  the  threshold  and  which  parts  are  larger. 
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^1  <  t  2 

•  points  with  intensity  <  t 

•  points  with  Intensity  >  1 1  and  <  t2 
o  points  with  Intensity  >  1 2 

(a) 


•  points  with  intensity  smaller  than  the  threshold 
o  points  with  intensity  larger  than  the  threshold 


(b) 


Figure  4.3:  (a)  Deriving  the  space  domain  constraint  for  the  two-dimensional  iterative  al¬ 
gorithm;  (b)  Ambiguous  range  information  via  samples  of  the  contour  associated  with  one 
threshold. 
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that  there  is  no  ambiguity  which  grid  points  are  inside  or  outside  of  a  given  threshold 


contour. 


2.  Let  jJ)  denote  the  value  of  /  in  the  /th  iteration  at  the  point  (^,  jft).  Then  the 

initial  guess  becomes: 

/( =  ±JuiW2)  B1jBj  =  0>  M  -  1 

M  M  2 

Similar  to  the  iterative  algorithm  of  section  (3.3.1),  the  convergence  of  this  algorithm  is 
independent  of  the  initial  guess. 

3.  Take  the  discrete  Fourier  transform  (DFT)  of  /  in  the  /th  step: 

ft'Hkukt)  =  0  *  <  M 


4.  Impose  the  bandlimited  constraint: 


FW(kltk2)  ,  0  <  *,,*2  <N  ,  M-N  <  kltk,  <  M 

0  ,  elsewhere 


5.  Inverse  discrete  Fourier  transform  of 


/(,+1)( =  DFT~1  [  ^(,+1)(fcl-fc2) 
M  M 


6.  Impose  the  space  domain  constraint: 

/(,+1)(T7.S)  =  *.(««,«.)  ,  /*,+l)(»,»)  <  *»(»!.»*) 


tu(ni»"2)  ,  /(i  +  1)(M’M)  >  *«(nl.na) 


7.  If  all  the  nodes  of  the  M  x  M  grid  satisfy  the  space  domain  constraint,  we  are  done. 
Otherwise  go  to  step  (3)  and  repeat  steps  (3)  through  (7). 


y 


The  above  algorithm  is  very  similar  to  the  one-dimensional  iterative  algorithm,  described 


in  section  (3.3.1).  Similar  to  the  one-dimensional  case,  ws  can  show  that  the  mean  square  error 
between  the  unknown  signal,  which  satisfies  the  space  and  frequency  domain  constraints,  and 


the  solution  obtained  via  successive  iterations  decreases.  Furthermore,  we  can  apply  Theorem 
(3.1)  to  prove  that  the  algorithm  is  guaranteed  to  converge  to  a  solution  which  satisfies  both 
the  space  and  frequency  domain  constraint,  if  such  a  solution  exists.  More  specifically,  the 
algorithm  described  in  this  section  is  a  special  case  of  the  iteration  described  in  Theorem  (3.1) 
where: 


The  Hilbert  space  H  is  the  finite  dimensional  space  of  all  M  x  M  point  real  sequences. 
The  inner  product  between  two  M  x  M  point  sequences  x(ni,n2)  and  y(ni,  n2)  is  given 


ij$* 

I 

vW, 


M- 1  M- 1 

(*,y)  =  Y  Y  x(ni>n2)  y(ni’n2) 

ni=0  ri2=0 

Ci  is  a  convex,  closed  subset  of  H,  and  includes  all  M  x  M  point  real,  bandlimited 
sequences  whose  DFT  is  0  in  the  range  N  <  k  <  M  -  N.  Furthermore,  C2  is  a  convex, 
closed  subset  of  H  containing  all  M  x  M  point  real  sequences,  y(ni,n2),  which  satisfy 
the  space  domain  constraint,  derived  from  the  level  crossing  information  as  described  in 
the  first  step  of  the  algorithm.  In  other  words,  any  M  x  M  sequence  which  satisfies 


m 

■&< 

m 


t«(ni.n2)  ^  y("i,n2)  <  tu(n1;n2),  V  0  <  ni,n2  <  Af 


belongs  to  the  set  C2.  Note  that  for  C2  to  be  closed,  it  is  necessary  to  define  it  exactly 
the  way  is  is  shown  above.  For  instance,  if  in  inequality  (4.8),  the  <  were  replaced  with 


1 


<,  then  C2  would  not  be  closed. 


Pi  is  the  bandlimiting  operator  which  operates  onto  the  set  C i,  and  P2  is  the  projection 
operator  which  imposes  the  space  domain  constraint  by  projecting  onto  C 2. 


•  The  relaxation  parameters  are  chosen  to  be  1. 


Thus,  the  iterative  algorithm  described  in  this  section,  is  a  special  case  of  the  iterative 
scheme  described  by  Theorem  (3.1).  Consequently,  it  converges  strongly  to  an  element  of 
Co  =  C\  n  C2,  the  set  of  two-dimensional  M  x  M  sequences  satisfying  both  the  space  and 
frequency  domain  constraints.  Since  in  practice  M  is  finite,  there  might  be  many  MxM  point 
real  sequences  in  the  set  Co-  However,  in  the  limit  as  M  — ♦  00,  our  theoretical  results  of  section 
(4.1)  guarantee  that  the  set  Co  will  contain  exactly  one  sequence.  Under  these  circumstances, 
H  will  not  be  finite  dimensional,  and  one  might  begin  to  wonder  whether  strong  convergence 
is  guaranteed  any  longer.  However,  as  it  is  mentioned  in  Theorem  (3.1),  since  at  least  one  of 
our  convex  closed  subsets,  namely  C 1,  is  finite  dimensional,  we  are  still  guaranteed  of  strong 
convergence  as  long  as  we  set  A 1  =  1. 


In  order  to  generalize  our  algorithm  to  the  case  where  the  relaxation  parameters  are  not 
equal  to  one,  we  must  modify  its  6th  step  in  the  following  way: 

•  Having  found  the  projection  of  ,  xf )  onto  Cj  given  by: 


/-■><2.=) -«!/«<=). g) 


apply  the  operator 


to  to  set; 


T\  —  1  +  Aj  (Pj  -  1) 


*<3 -5)  - '' "<3 -S,  +  *•  1  ’(S-g) 
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•  Now  impose  the  space  domain  constraint  by  projecting  /^+1KlJ>  m )  onto  ^2  to  8et: 

9.9)  a  al/“+,)(S».9)l 

/*,+1)(».S)  ,  <  /,,+1,(9.&)  <<«(»,,»,) 

=  M-l.l)  ,  /"+,»(9,»)  < 

Mni>ns)  .  /(,+I)(9>9)  >  *u(ni."i) 

•  Apply  the  operator  T2  given  by: 

T2  =  1  +  A2  (  P2  -  1) 

‘«/l'+1|(9.9)ws«: 

/(,+,|(9.9)  =  rI|/(,t,,(9.9>l 

=  /|,+l)(».»)  +  W,+,,(9.9>  -  /t,+l,(9.»)) 

The  main  advantage  of  the  generalized  version  of  the  iterative  algorithm  is  that  its  convergence 
properties  can  be  modified  by  adjusting  the  relaxation  parameters  Aj  and  A2.  In  fact,  there 
has  been  great  deal  of  research  on  the  “optimum”  choice  of  relaxation  parameters;  we  will  only 
mention  the  simplest  such  technique,  over-relaxation,  which  essentially  consists  of  choosing  the 
relaxation  parameters  in  the  range  (1,2).  We  have  found  experimentally  that  convergence  rates 
do  not  present  serious  enough  problems  to  seek  other  means  of  accelerating  convergence. 

4.2.2  Examples  of  the  Iterative  Algorithm 

An  example  of  the  reconstruction  of  the  original  eye  picture  of  figure  (2.10),  is  shown  in 
figures  (4.4)  and  (4.5).  Figure  (4.4)  shows  the  5  level  non-uniform  amplitude  quantized  version 
of  figure  (2.10)  and  the  threshold  contours  associated  with  these  levels.  The  reconstructed 
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version  of  figure  (4.4)  is  shown  in  figure  (4.5);  the  quantizing  grid  was  256  x  256,  relaxation 
parameters  were  Ai  =  1  and  A2  =  1.7,  and  the  number  of  iterations  was  20.  We  have  found 
experimentally  that  in  order  to  reconstruct  the  image  from  fewer  than  5  thresholds,  the  grid 
size  has  to  be  much  finer  than  512  x  512.  We  will  investigate  the  trade  off  between  number 
of  thresholds  and  the  size  of  the  quantizing  grid  in  more  detail  in  Chapter  5.  In  this  section, 
however,  we  will  only  show  three  examples  of  reconstruction  for  different  number  of  thresholds 
and  various  values  of  M.  Figures  (4.6),  (4.7),  and  (4.8)  show  the  reconstructed  version  of 
the  eye,  cman  and  vegas  picture  from  4,  6  and  8  thresholds  via  the  iterative  algorithm  for 
different  values  of  M.  The  relaxation  parameters  for  all  the  reconstructions  were  chosen  to 
be  Aj  =  1,  A2  =  1.75.  As  shown,  increasing  the  number  of  thresholds  affects  the  quality  of 
reconstruction  in  a  more  substantial  way  than  increasing  M.  In  addition,  comparing  figures 
(4.6)  and  (3.17)  we  can  see  that  the  quality  of  reconstruction  via  the  two-dimensional  iterative 
algorithm  of  this  section  is  superior  to  the  one-dimensional  iterative  algorithm  of  section  (3.3.1). 


Similar  to  the  one-dimensional  case,  we  have  found  experimentally  that  altering  the  6th 
step  of  the  algorithm  in  the  following  manner: 


/(i+1)(^,^)  .  <  /(<+1)(M)  <  *«(»!,»») 


Mm.n?)  +  U(ni,na) 
2 


otherwise 


reduces  the  mean  square  error  faster  than  the  generalized  version  of  the  iterative  algorithm 
where  relaxation  parameters  are  used  to  accelerate  convergence.  This  is  true,  in  spite  of  the 
fact  that  the  mean  square  error  does  not  necessarily  decrease  from  one  iteration  to  the  next  and 
that  the  convergence  of  this  technique  can  not  even  be  guaranteed  by  applying  Theorem  (3.1). 
Although  this  faster  convergence  can  be  explained  by  formulating  the  problem  in  a  statistical 
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framework  [49],  we  will  not  develop  it  here.  The  important  point  is  that,  although  in  most  re¬ 
constructions  carried  out  via  the  above  modified  version,  the  algorithm  converges  to  a  sequence 
satisfying  both  the  space  and  frequency  domain  constraints,  the  quality  of  the  reconstructed 
image  was  almost  identical  to  the  one  obtained  via  the  over  relaxation  method.  This  suggests 
that  the  reason  behind  poor  quality  of  reconstruction  for  small  number  of  thresholds  is  not  the 
fact  that  the  algorithm  did  not  converge,  but  that,  as  with  the  semi-implicit  iterative  scheme, 
there  are  many  signals  satisfying  the  space  and  frequency  domain  constraints,  and  the  algo¬ 
rithm  has  simply  converged  to  one  of  them.  Increasing  the  number  of  thresholds  or  increasing 
M  imposes  more  constraints  on  the  signal,  thereby  reducing  the  number  of  sequences  satisfying 
both  the  space  and  frequency  domain  requirements.  As  we  said  earlier,  the  solution  via  the 
iterative  algorithm  becomes  unique  as  M  — *  oo. 

Before  we  end  this  section,  let  us  briefly  compare  the  iterative  and  linear  least-squares  ap¬ 
proach  for  reconstruction  via  the  implicit  sampling  strategy,  described  in  this  chapter.  Unlike 
the  iterative  approach,  which  needs  all  the  quantized  threshold  contours,  the  linear  least-squares 
approach  only  needs  samples  of  the  threshold  contours.  While,  the  iterative  algorithm  is  less 
storage  and  computation  intensive,  the  quality  of  reconstruction  via  the  iterative  algorithm  is 
not  as  high  as  the  least-squares  approach.  However,  one  must  keep  in  mind  the  fact  that  recon¬ 
structions  using  the  latter  technique  use,  at  least  in  principle,  the  exact  location  of  the  crossings, 
whereas  the  iterative  algorithm  uses  the  rather  coarsely  quantized  version  of  the  crossings,  as 
long  as  M  is  finite.  Thus,  in  situations  where  the  location  of  the  level  crossings  is  quantized  too 
coarsely,  unlike  the  linear  least-squares  algorithm,  which  is  completely  unsuccessful  in  carrying 
out  reconstruction,  the  iterative  algorithm  is  capable  of  producing  approximate  reconstructions. 
We  will  study  these  quantization  issues  at  length  in  Chapter  5. 
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Figure  4.8:  Reconstructed  version  of  the  vegas  picture  via  the  two-dimensional  iterative  algo¬ 
rithm  of  this  section.  The  number  of  thresholds  are  4,6  and  8,  increasing  from  left  to  right, 
and  the  number  of  equally  spiked  samples  on  lines,  M  are  127  (no  reconstruction),  128,  256 
increasing  from  top  to  bottom. 


to 


Chapter  5 


Preliminary  Speculations  on 
Quantization  Properties 


In  the  past  three  chapters,  we  have  proposed  various  approaches  to  the  problem  of  recon¬ 
struction  from  level  crossings.  As  we  mentioned  in  Chapter  1,  our  main  motivation  for  solving 
this  problem  has  been  to  find  sampling  strategies  whose  characteristics  lie  in  between  the  ex¬ 
plicit  Nyquist  sampling  and  the  implicit  zero  crossing  strategy  as  defined  by  Curtis  [51].  More 
specifically,  in  Nyquist  sampling,  the  amplitude  of  the  signal  is  specified  to  infinite  precision 
at  prespecified  points,  and  all  the  bits  used  to  represent  the  signal  are  essentially  amplitude 
bits.  On  the  other  hand,  for  successful  reconstruction  from  zero  crossings,  the  position  of  the 
crossings  need  to  be  known  extremely  accurately,  whereas  only  one  bit  is  needed  to  quantize  the 
amplitude  information.  Similarly,  while  the  bandwidth  used  in  representing  signals  via  their 
Nyquist  samples  is  minimal  and  the  dynamic  range  is  maximal,  the  zero  crossing  representation 
requires  maximal  bandwidth  and  minimal  dynamic  range.  As  will  see  in  this  chapter,  in  repre¬ 
senting  signals  via  their  samples  of  multiple  level  threshold  crossings,  the  required  bandwidth 
and  dynamic  range  are  in  between  those  of  the  zero  crossing  and  Nyquist  representations. 

To  demonstrate  this  more  quantitatively,  we  need  to  investigate  position  and  amplitude 

quantization  requirements  of  our  sampling/reconstruction  schemes  as  a  function  of  the  number 
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of  thresholds.  A  rigorous  and  thorough  investigation  of  these  quantization  characterizations 
involves  extensive  dealing  with  coding  issues  and  experiments,  both  of  which  are  beyond  the 
scope  of  this  thesis.  Therefore,  the  nature  of  discussions  in  this  chapter  tends  to  be  rather 
preliminary,  and  most  of  the  conclusions  are  somewhat  tentative.  However,  these  speculative 
results  can  be  used  as  a  starting  point  for  further  research  in  the  areas  of  multidimensional 
signal  representation  and  image  coding. 

Since  we  have  proposed  numerous  sampling  and  reconstruction  strategies  in  the  last  three 
chapters,  even  a  preliminary  investigation  of  their  quantization  properties  is  a  formidable  task. 
Thus,  our  approach  in  this  chapter  has  been  to  study  a  subset  of  these  schemes.  More  specifi¬ 
cally,  we  will  study  the  quantization  properties  of  the  linear  least-squares  approach  in  section 
(5.1),  and  the  iterative  algorithms  in  section  (5.2).  The  reason  behind  choosing  the  linear  least- 
squares  approach  is  that,  if  it  is  implemented  via  stable  algorithms  such  as  QR  decomposition,  it 
can  result  in  extremely  robust  reconstructions  for  both  the  semi-implicit  and  implicit  sampling 
strategies.  The  main  reason  for  studying  the  iterative  algorithms  is  that  their  quantization 
characteristics  are  substantially  different  from  our  other  algorithms. 

The  organization  of  this  chapter  is  as  follows.  In  sections  (5.1)  and  (5.2),  we  will  examine 
the  quantization  requirements  of  the  linear  least-squares  approach  and  iterative  algorithms  as 
a  function  of  the  number  of  thresholds  and  the  sampling  strategy.  As  we  will  see  in  section 
(5.1),  the  number  and  geometric  distribution  of  the  level  crossings  play  an  important  role  in 
the  robustness  of  linear  least-squares  reconstruction  via  the  QR  decomposition.  In  addition, 
for  both  the  QR  decomposition  and  iterative  algorithms,  increasing  the  number  of  thresholds, 
which  in  effect  corresponds  to  an  increase  in  the  number  of  amplitude  bits,  leads  to  a  decrease  in 
the  required  number  of  position  bits.  Thus,  there  is  a  tradeoff  between  the  number  of  thresholds 


and  total  number  of  amplitude  and  position  bits.  Finally,  in  section  (5.3),  we  show  that  under 
certain  circumstances  Nyquist  sampling  simply  becomes  a  special  case  of  our  semi-implicit  and 
implicit  sampling  strategies;  this  will  bridge  the  conceptual  gap  between  explicit,  semi-implicit 
and  implicit  sampling  strategies,  unify  seemingly  unrelated  sampling  schemes  and  provide  us 
with  a  spectrum  of  sampling  techniques  for  multidimensional  signals. 

5.1  Linear  Least-Squares  Approach 

This  section  includes  a  preliminary  investigation  of  the  quantization  properties  of  recon¬ 
struction  via  the  linear  least-squares  approach  and  the  QR  decomposition.  Before  approaching 
this  problem  however,  we  need  to  address  a  few  issues. 

The  first  issue  has  to  do  with  the  fact  that  quantization  procedures  for  the  implicit  and 
semi-implicit  sampling  approaches  are  somewhat  different  from  each  other.  As  we  will  see  in 
section  (5.1.1),  while  for  the  implicit  approach,  level  crossing  samples  are  quantized  on  a  square 
grid,  for  the  semi-implicit  approach,  we  exploit  the  geometric  constraint  of  the  sample  positions 
to  quantize  their  locations. 

The  second  issue  that  needs  to  be  addressed  is  how  the  number  of  reconstruction  samples 
influences  reconstruction  robustness.  As  we  will  see  in  section  (5.1.2),  increasing  the  number 
of  samples  will  result  in  fewer  position  bits  per  sample,  although  beyond  a  certain  level,  it 
increases  the  total  number  of  bits  used  to  represent  a  signal. 

Having  discussed  these  two  issues  in  sections  (5.1.1)  and  (5.1.2),  in  section  (5.1.3),  we  will 
investigate  how  the  number  of  thresholds  affect  position  and  amplitude  quantization  charac 
teristics  of  reconstruction  via  the  linear  least-squares  approach  and  the  QR  decomposition. 
Our  experimental  results  seem  to  suggest  that,  although  increasing  the  number  of  thresholds 
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Figure  5.1:  Quantization  procedure  for  the  implicit  sampling 
strategy  via  the  linear  least-squares  approach. 


initially  decreases  the  number  of  position  bits,  beyond  a  certain  point,  it  increases  the  total 
number  of  position  and  amplitude  bits  required  for  representing  a  two  dimensional  signal. 
5.1.1  Quantization  Procedures 

Quantization  procedures  which  we  used  for  the  implicit  and  semi-implicit  sampling  strategies 
are  slightly  different  from  each  other.  Figure  (5.1)  shows  quantization  of  samples  of  threshold 
contours  obtained  via  the  implicit  sampling  approach.  As  shown,  we  used  a  simple  quantization 
algorithm  in  which  the  quantized  coordinates  of  a  given  sample  are  chosen  to  be  the  coordinates 
of  the  center  point  of  the  square  the  sample  falls  in.  In  addition,  for  situations  in  which  two 
or  more  samples  fall  into  the  same  quantizing  square,  the  sample  closest  to  the  center  of  the 


square  is  kept  and  the  remaining  ones  are  discarded. 


For  the  semi-implicit  approach,  we  can  take  advantage  of  the  geometry  of  the  sampling  lines. 
More  specifically,  in  a  sampling  scenario  with  nj  sampling  lines,  our  strategy  for  specifying  the 
location  of  a  given  sample  has  been  to: 

1.  Use  log2  ni  bits  to  specify  the  line  it  falls  on. 

2.  Use  b  bits  to  specify  the  location  of  sample  along  its  sampling  line. 

It  is  worthwhile  to  point  out  that  there  are  many  other  ways  to  quantize  the  position  of  our 
samples,  and  our  proposed  coding  schemes  are  almost  certainly  not  optimal.  In  fact,  a  rigorous 
investigation  of  the  quantization  properties  involves  dealing  with  coding  issues  which  are  well 
beyond  the  scope  of  this  thesis.  Since  we  have  chosen  unsophisticated  and  simple  quantization 
strategies,  it  is  important  to  bear  in  mind  that  the  nature  of  our  conclusions  tend  to  be  rather 
speculative  and  preliminary. 

5.1.2  Choice  of  the  Number  of  Samples 

Although  the  theoretical  results  of  Chapters  2  and  4,  provide  us  with  a  variety  of  sufficient 
conditions  under  which  samples  of  level  crossings  of  a  BLP  signal  can  be  used  to  uniquely 
specify  it,  using  more  interpolation  points  than  the  minimum  required  by  these  theoretical 
results  should  improve  the  robustness  of  reconstruction,  and  thus  decrease  the  required  number 
of  position  bits  for  specifying  the  location  of- the  crossings.  Although  for  a  fixed  quality  of 
reconstruction,  increasing  the  number  of  samples  will  initially  decrease  the  number  of  position 
bits  per  sample,  beyond  a  certain  point,  it  increases  the  total  number  of  bits  used  to  represent 
a  signal.  Thus,  our  goal  in  this  section  is  to  find  the  extent  to  which  oversampling  signals  will 
reduce  the  total  number  of  bits  necessary  for  representing  them.  To  this  end,  we  carried  out 
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Figure  5.2:  The  original  eye.lp  picture  which  is  the  low  pass 
version  of  the  31  x  31  eye  picture. 

a  series  of  experiments  on  a  picture  shown  in  figure  (5.2)  with  15  x  15  region  of  support  in 
the  Fourier  domain.  This  picture  which  is  more  or  less  the  low  pass  version  of  the  original  eye 
picture  of  figure  (2.10),  will  be  referred  to  as  eye.lp. 

For  the  semi-implicit  approach,  we  first  found  all  the  intersections  of  7  level  crossing  contours 
with  15  equidistant  lines  of  unity  slope,  and  then  chose  five  subsets  of  these  intersections, 
with  a  different  number  of  samples  in  each  subset.  To  guarantee  unique  reconstructibility, 
the  distribution  of  the  interpolation  points  among  various  sampling  lines  for  each  subset  was 
chosen  in  such  a  way  that  the  conditions  of  Theorem  (2.3)  were  satisfied.  As  it  turns  out, 
for  a  small  number  of  thresholds,  the  quality  of  recovered  images  via  the  QR  decomposition 
deteriorates  very  rapidly  once  the  number  of  bits  used  for  quantizing  the  location  of  samples 

drops  below  a  certain  level.  Hence,  it  is  rather  straightforward  to  find  the  minimum  number  of 
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position  bits  required  for  successful  reconstruction.  Nevertheless,  it  is  important  to  emphasize 
that  our  evaluation  of  the  reconstructed  images  were  purely  subjective  and  rather  informal.  A 
plot  of  the  minimum  number  of  position  bits  per  sample  versus  the  number  of  reconstruction 
samples  is  shown  in  figure  (5.3).  As  shown,  the  number  of  position  bits  drops  very  rapidly  as 
the  number  of  samples  is  increased  from  225  to  250.  However,  this  drop  becomes  considerably 
less  significant  as  the  number  of  samples  is  further  increased  to  399.  Figure  (5.3)  also  contains 
the  plot  of  the  position  bits  versus  the  number  of  samples  for  the  implicit  sampling  approach  of 
Chapter  4.  Like  the  semi-implicit  approach,  the  number  of  position  bits  drops  considerably  as 
the  number  of  samples  is  increased  from  225  to  256;  however,  increasing  the  number  of  samples 
from  256  to  576  does  not  seem  to  decrease  the  number  of  position  bits  any  further. 

The  shape  of  the  curves  in  figure  (5.3)  suggests  that  a  plot  of  total  number  of  position 
bits  versus  the  number  of  samples  must  exhibit  a  minimum.  Indeed,  as  shown  in  figure  (5.4), 
the  minimum  occurs  when  the  number  of  samples  is  around  250.  It  is  important  to  note  that 
the  vertical  axis  in  figure  (5.3)  indicates  the  number  of  position  bits  per  sample,  whereas  in 
figure  (5.4)  it  stands  for  the  total  number  of  position  bits  normalized  to  the  number  of  Fourier 
coefficient  i.e.: 

number  of  position  bits  per  sample  x  number  of  samples 
number  of  Fourier  coef  ficients 

The  curves  in  figure  (5.4)  suggest  that  for  a  fixed  number  of  thresholds,  (in  this  case  7), 
oversampling  the  level  crossings  by  approximately  10  percent  results  in  lowest  total  number  of 
position  bits.  In  the  next  section,  we  will  investigate  how  the  number  of  thresholds  and  the 
number  of  reconstruction  samples  affect  the  quantization  requirements  of  reconstruction  via 
the  QR  decomposition. 
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Figure  5.3:  Number  of  position  bits  per  sample  used  for  successful  reconstruction  of  the  eye. Ip 
picture  versus  the  number  of  samples.  Each  curve  represents  more  or  less  constant  quality  as 
judged  subjectively  and  informally. 
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Figure  5.4:  Normalized  number  of  position  bits  used  for  successful  reconstruction  of  the  eye. Ip 
picture  versus  the  number  of  samples.  Each  curve  represents  more  or  less  constant  quality  as 
judged  subjectively  and  informally. 


5.1.3  Quantization  Requirements 

As  we  mentioned  earlier,  position  and  amplitude  quantization  requirements  of  various  re¬ 
construction  strategies  depend  on  factors  such  as  sampling  strategy  and  the  number  of  recon¬ 
struction  samples.  Having  discussed  the  effect  of  number  of  reconstruction  samples  in  section 
(5.1.2),  we  are  now  ready  to  investigate  quantization  requirements  of  semi-implicit  and  implicit 
sampling  strategies  as  a  function  of  the  number  of  thresholds. 

Figure  (5.5)  shows  the  plot  of  mean  square  error  (mse)  versus  the  total  number  of  amplitude 
and  position  bits  as  a  function  of  the  number  of  thresholds.  The  mean  square  error  between 
the  original  and  reconstructed  image  is  given  by: 

mae  =  iiwTTP  £  £  inkuk,)  -  F(fc 

V  '  fc,  =  -N  ks  =  -N 

where  F(k\,  fc2)  and  F(ki,  k2)  correspond  to  the  Fourier  coefficients  of  the  original  and  recon¬ 
structed  image  respectively.  We  have  found  experimentally  that  the  quality  of  the  reconstructed 
image  becomes  almost  indistinguishable  from  the  original  one  when  mse  <  .1  .  Figure  (5.5b) 
corresponds  to  the  implicit  sampling  strategy  and  (5.5a)  corresponds  to  semi-implicit  sampling 
strategy  with  equidistant  lines  of  unity  slope.  In  both  cases,  the  number  of  samples  used  for 
reconstruction  was  225  i.e.  the  minimum  number  required  by  our  theoretical  results.  As  we 
would  expect,  the  slope  of  the  curves  shown  in  figure  (5.5)  are  negative,  indicating  that  the 
quality  of  the  reconstruction  is  improved  as  the  the  number  of  position  bits  is  increased.  In 
addition,  the  spacing  between  the  curves  decreases  from  right  to  left,  indicating  that  the  im¬ 
provement  in  the  quality  of  reconstruction  decreases  as  the  the  number  of  thresholds  gets  larger. 
In  fact,  as  shown  in  figure  (5.5b),  for  mse  <  .05  the  quantization  characteristics  of  the  curve 
corresponding  to  16  thresholds  is  slightly  more  favorable  than  that  of  32  thresholds. 


Superimposing  curves  of  figures  (5.5a)  and  (5.5b),  we  obtain  figure  (5.6)  which  shows  the 
semi-implicit  and  implicit  curves  corresponding  to  7,  16,  and  32  thresholds.  As  shown,  the  slope 
of  the  semi-implicit  curves  are  more  negative  than  the  slope  of  implicit  ones.  In  addition,  it  ap¬ 
pears  that  for  mse  >1,  the  implicit  curves  exhibit  more  favorable  quantization  characteristics 
than  the  semi-implicit  ones.  This  is  primarily  due  to  the  fact  that  for  semi-implicit  sampling 
situations  in  which  the  number  of  reconstruction  samples  does  not  exceed  the  required  mini¬ 
mum  value,  the  geometric  distribution  of  samples  is  constrained  by  the  conditions  of  theorems 
of  Chapter  2,  whereas  the  interpolation  points  of  the  implicit  sampling  are  more  or  less  chosen 
at  random  and  can  be  more  evenly  distributed  across  the  image.  As  an  example,  consider 
semi-implicit  sampling  scenarios  where  all  the  sampling  lines  have  identical  integer  slopes  m. 
In  this  case,  the  ratio  between  the  maximum  and  minimum  number  of  required  samples  on  the 
sampling  lines  is  given  by: 

__  maximum  number  of  points  on  a  line  _  1  +  2(|m|  +  1  )  N  (5  1) 

^  minimum  number  of  samples  on  a  line  1  +  2(|m|  +  1)  N  -  4N  K  ‘  ' 


which  can  be  approximated  by: 


H  +  1 

Iml  -  1 


Clearly,  an  even  distribution  of  samples  corresponds  to  small  values  of  p.  Therefore,  for  |m|  >  1 
increasing  |m|,  and  for  |m|  <  1,  decreasing  |m[  will  result  in  smaller  value  of  p,  and  thus,  a 
more  favorable  distribution  of  samples.  To  verify  this  experimentally,  a  reconstructed  version 
of  the  eye. Ip  picture  using  samples  on  lines  of  slopes  3  and  1  are  shown  in  figure  (5.7).  The 
artifacts  in  figure  (5.7a)  occur  at  a  place  where  two  neighboring  sampling  lines  of  unity  slope 
have  15  and  1  points  respectively.  Indeed,  the  value  of  p  is  15  for  sampling  lines  of  unity  slope 
and  approximately  2  for  lines  of  slope  3. 
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Figure  5.6:  Plot  of  mse  versus  total  normalized  number  of  position  and  amplitude  bits  for 
semi-implicit  and  implicit  sampling.  The  number  of  reconstruction  samples  was  fixed  at  225 
and  the  number  of  thresholds  was  7  in  (a),  16  in  (b),  and  32  in  (c). 


As  the  number  of  reconstruction  samples  is  increased  beyond  the  minimum  required  by  our 
theoretical  results,  the  distribution  of  interpolation  points  becomes  more  even,  thus  resulting 
in  fewer  position  bits.  Figure  (5.8)  shows  the  plot  of  mse  versus  total  number  of  position  and 
amplitude  bits  as  a  function  of  the  number  of  thresholds.  The  curves  shown  in  figure  (5.8)  were 
generated  in  a  similar  fashion  to  those  of  figure  (5.5)  except  that  the  number  of  reconstruction 
samples  was  256  instead  of  225.  Notice  that  for  the  semi-implicit  curves  of  figure  (5.8a), 
increasing  the  number  of  thresholds  from  7  to  16  improves  the  quantization  characteristics, 
while  further  increase  from  32  to  64  thresholds  degrades  it.  As  far  as  implicit  sampling  curves 
of  figure  (5.8b)  go,  the  “optimum”  number  of  thresholds  which  results  in  minimum  number  of 
bits  varies  as  a  function  of  mse.  For  instance,  for  .2  <  mse  <  1  it  is  7,  for  .06  <  mse  <  .2  it  is 
16,  and  for  .02  <  mse  <  .06  it  is  32. 

Although  it  is  rather  difficult  to  draw  strong  conclusions  from  figures  (5.5)  and  (5.8),  the 
plots  seem  to  suggest  that  for  both  semi-implicit  and  implicit  sampling  strategies,  there  exists 
an  “optimum”  number  of  thresholds  which  results  in  the  least  number  of  amplitude  and  position 
bits.  This  “optimum”  number  appears  to  be  a  function  of  mse,  the  sampling  strategy  and  the 
number  of  reconstruction  samples. 

5.2  Iterative  Approach 

As  we  mentioned  earlier,  the  quantization  properties  of  the  iterative  algorithms  of  sections 
(3.3.1)  and  (4.2.1)  are  substantially  different  from  our  other  algorithms.  In  this  section,  we  will 
present  some  preliminary  results  on  the  quantization  properties  of  these  iterative  algorithms. 
More  specifically,  sections  (5.2.1)  and  (5.2.2)  include  amplitude  and  position  quantization  re- 


quirements  of  the  semi-implicit  and  the  implicit  sampling  strategies  respectively. 

5.2.1  Quantization  Requirements  for  Semi-implicit  Sampling 

In  this  section,  we  will  present  a  preliminary  investigation  of  the  quantization  properties  of 
the  iterative  algorithm  for  the  semi-implicit  sampling  strategy.  Figure  (5.9)  shows  the  mean 
square  error  between  the  original  eye. Ip  picture  and  its  reconstructed  version,  versus  the  number 
of  position  and  amplitude  bits.  The  number  of  thresholds  associated  with  the  four  curves  shown 
in  figure  (5.9)  is  6,  8,  12,  and  16.  For  each  curve  the  thresholds  were  chosen  with  equal  spacing 
between  0  and  256.  The  four  points  on  each  curve  correspond  to  different  number  of  equally 
spaced  points  on  the  sampling  lines,  i.e.  M  —  32,  64,  128.  In  addition,  the  sampling  lines 
were  chosen  to  be  equidistant  and  of  unity  slope.  Although  our  theoretical  results  of  Chapter  2 
only  require  15  sampling  lines  of  unity  slope,  to  be  able  to  use  the  FFT  with  a  power  of  2  for  the 
interpolation  part  of  the  iterative  algorithm  (i.e.  interpolation  from  recovered  one-dimensional 
signals  to  a  square  grid  on  the  two-dimensional  one),  we  chose  16  equally  spaced  sampling 
lines.  The  y  axis  in  figure  (5.9)  corresponds  to  mse,  and  for  completeness  we  have  also  included 
the  actual  reconstructed  images  for  different  values  of  M  and  different  number  of  thresholds 
in  figure  (5.10).  As  shown  in  figure  (5.10),  the  quality  among  the  reconstructed  images  with 
approximately  equal  mse  is  comparable.  Also,  the  quality  of  reconstruction  certainly  improves 
as  mse  is  decreased.  Finally,  the  x  axis  in  figure  (5.9)  indicates  the  total  number  of  position  and 
amplitude  bits  used  for  the  particular  reconstruction  at  hand.  Because  of  the  inherent  structure 
of  the  samples  used  by  the  iterative  algorithm,  there  are  a  variety  of  ways  to  represent  the  signal 
under  consideration  and  to  arrive  at  the  number  of  quantization  bits.  The  most  straightforward 
way  is  to  quantize  the  location  of  threshold  crossings  on  the  sampling  lines  to  log2  M  bits,  so 
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5.9:  Plot  of  mean  square  error  versus  normalized  number  of  position  and  amplitude 
ation  bits  as  a  function  of  the  number  of  thresholds  for  the  eye. Ip  picture. 
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that  the  space  domain  constraint  for  the  M  equally  spaced  points  on  each  of  the  sampling  lines 
can  be  derived  easily.  Using  this  strategy,  the  total  number  of  amplitude  and  position  bits  for 
representing  the  signal  can  be  written  as: 


(log2  M  +  log2nt  +  log2ni)  ^  IV,- 

*=i 


where 


•  Nt  denotes  the  number  of  intersections  of  the  tth  threshold  crossings  with  all  the  sampling 
lines. 


•  n(  denotes  the  number  of  thresholds. 

•  r»i  denotes  the  number  of  sampling  lines. 

•  M  denotes  the  number  of  equally  spaced  points  on  each  of  the  sampling  lines. 


An  alternative  strategy  for  representing  the  signal  is  to  specify  the  range  of  intensity  for  each  of 
the  equally  spaced  points  on  the  sampling  lines.  More  specifically,  if  the  number  of  thresholds 
is  nt,  and  the  signal  is  known  to  be  in  the  range  [0,256],  then  the  value  of  the  signal  at  any 
given  point  lies  in  one  of  the  (nt  +  1)  intervals  corresponding  to  our  nt  thresholds.  In  this  case, 
the  total  number  of  bits  used  to  represent  the  signal  is  given  by  Mnj  log2(nt  +  1).  Clearly,  this 
second  quantization  strategy  outperforms  the  first  one  for  large  values  of  and  small  values  of 
M .  Our  strategy  in  determining  the  total  number  of  position  bits  for  the  abscissa  of  figure  (5.9) 
has  been  to  choose  the  minimum  of  the  above  two  quantization  strategies.  As  it  turns  out,  we 
would  reach  the  same  conclusions  if  we  choose  either  of  the  quantization  strategies  described 
above. 


Having  discussed  the  details  of  generating  the  curves  shown  in  figure  (5.9)  and  the  images 
of  figure  (5.10),  it  is  now  appropriate  to  make  a  few  observations  and  comments  regarding 


their  shapes.  As  we  would  expect,  the  slope  of  each  curve  is  negative  indicating  that  for  fixed 
number  of  thresholds  the  mean  square  error  decreases  as  M  is  increased.  In  addition,  for  fixed 
M,  the  mean  square  error  is  decreased  as  the  number  of  thresholds  is  increased.  The  decreasing 
distance  between  the  curves  shown  in  figure  (5.9)  is  indicative  of  the  fact  that  as  the  number 
of  thresholds  increases,  the  resulting  drop  in  mse  is  less.  Thus,  the  decrease  in  mse  as  the 
number  of  thresholds  changes  from  6  to  8  is  more  substantial  than  when  it  changes  from  12 
to  16.  An  interesting  question  to  address,  however,  is  whether  or  not  there  is  an  “optimum” 
number  of  thresholds  for  which  the  lowest  number  of  amplitude  and  position  quantization  bits 
is  achieved.  As  figure  (5.9)  shows,  this  “optimum”  number  varies  as  a  function  of  the  mean 
square  error.  For  instance,  for  the  value  of  mse  in  the  range  [0.53, 0.85],  it  is  between  6  to  8,  and 
for  mse  in  the  range  [0.15,0.23],  it  is  between  12  and  16.  Thus,  for  iterative  reconstruction  via 
semi-implicit  sampling,  it  appears  that  the  “optimum”  number  of  thresholds,  resulting  in  least 
number  of  amplitude  and  position  bits,  is  a  decreasing  function  of  mse.  In  the  next  section,  we 
will  investigate  iterative  reconstruction  via  implicit  sampling. 

5.2.2  Quantization  Requirements  for  Implicit  Sampling 

Our  main  goal  in  this  section  is  to  investigate  quantization  properties  of  the  iterative  recon¬ 
struction  algorithm  for  the  implicit  sampling  strategy.  The  basic  idea  behind  this  algorithm,  is 
to  iterate  between  the  space  domain  constraints  derived  from  the  quantized  threshold  contours, 
and  the  frequency  domain  constraints  derived  from  the  bandlimitedness  of  the  signal. 

Similar  to  the  iterative  reconstruction  algorithm  of  the  semi-implicit  sampling,  there  are 
several  strategies  one  might  take  to  arrive  at  the  total  number  of  amplitude  and  position 
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quantization  b:tc  for  representing  a  given  image  via  the  implicit  sampling  scheme.  Since  in 
each  of  these  strategies,  we  are  encoding  the  quantized  contours  corresponding  to  different 
thresholds,  almost  all  of  the  efficient  coding  schemes  for  two  tone  images  proposed  and  studied 
by  many  researchers  [54,55]  can  be  used  for  representing  our  images.  The  most  obvious  way  of 
encoding  the  boundary  points  of  a  threshold  contour  quantized  to  6  bits  is  to  use  26  bits  for 
x  and  y  coordinates  of  each  point.  Using  this  strategy,  the  total  number  of  bits  required  to 
specify  n(  threshold  contours  is  given  by : 

(log2  nt  +  26)(^ATi(6)) 

«=i 

where  JV,(6)  denotes  the  number  of  quantized  boundary  points  on  a  26  x  2b  quantization  grid 
for  the  ith  threshold.  A  second  and  more  efficient  way  of  representing  the  boundary  points  is  to 
follow  the  boundary  i.e.  to  do  contour  tracing.  Since  the  image  is  quantized  on  a  2l  x  26  square 
grid,  each  pel  has  only  eight  neighbors.  Therefore,  it  is  sufficient  to  use  3  bits  to  indicate  where 
the  next  boundary  point  is.  Of  course,  to  get  on  each  boundary,  we  need  to  specify  the  position 
of  an  initial  point.  Thus,  ignoring  the  additional  bits  required  to  get  to  the  initial  points  on 
the  boundaries  and  to  specify  their  associated  threshold  value,  total  number  of  points  required 
for  specifying  an  image  is  given  by 

3  £>,(6) 

i=i 

Our  third  encoding  strategy  involves  specifying  the  range  of  the  signal  for  each  node  of  the 
2*  x  2*  quantization  grid.  This  strategy  is  suitable  for  sampling  schemes  with  large  number 
of  thresholds  and  coarse  quantization  of  the  contours.  The  total  number  of  bits  required  for 
specifying  an  image  via  this  encoding  technique  is  given  by 

226  log2(n(  4-  1) 
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As  far  as  efficiency  of  these  encoding  strategies  goes,  clearly  the  second  strategy  outperforms 
the  first  one,  and  the  relative  efficiencies  of  the  second  and  third  encoding  schemes  depend  on 
the  form  of  JVj( b).  Intuitively,  we  would  expect  Ni(b)  to  be  proportional  to  2b,  since  the  number 
of  boundary  points  of  quantized  threshold  contours  is  doubled  as  the  size  of  quantizing  grid  is 
increased  from  2*  x  26  to  2*+1  x  2b+l.  Thus,  for  small  values  of  b  and  large  values  of  nt,  the 
third  strategy  outperforms  the  second  one.  As  it  turns  out,  our  major  conclusions  are  more 
or  less  independent  of  the  actual  encoding  strategy  used;  so  our  adopted  strategy  has  been  to 
choose  the  minimum  of  second  and  third  quantization  strategies  for  representing  images. 

Having  described  our  quantization  strategy,  we  will  now  examine  how  the  number  of  thresh¬ 
olds  affect  the  required  number  of  position  bits  and  the  quality  of  reconstruction.  Figure  (5.11) 
shows  the  plot  of  the  mean  square  error  between  the  eye. Ip  picture  and  its  reconstructed  version 
versus  the  total  number  of  position  and  amplitude  bits.  The  four  curves  shown  in  figure  (5.11) 
correspond  to  four  different  values  of  the  grid  size,  i.e  16, 32,64  and  128.  Various  points  on  each 
curve  correspond  to  reconstruction  from  different  numbers  of  thresholds.  The  reconstructed  im¬ 
ages  corresponding  to  different  points  on  the  curves  of  figure  (5.11)  are  shown  in  figure  (5.12). 
The  reconstruction  was  carried  out  via  the  iterative  algorithm  of  section  (4.2)  with  Aj  =  1.5 
and  A2  =  1.95.  These  values  of  A’s  were  found  to  result  in  the  smallest  values  of  the  mean 
square  error.  The  following  observations  can  be  made  about  the  plots  shown  in  figure  (5.11): 

•  The  slopes  of  the  curves  are  negative  indicating  that  the  quality  of  reconstruction  is 
improved  as  the  number  of  thresholds  is  increased  or  equivalently  as  the  number  of  position 
and  amplitude  bits  used  for  representing  the  image  is  increased. 

•  Smallest  number  of  position  and  amplitude  bits  is  achieved  when  the  grid  size  is  at  its 
minimum  value  i.e.  16,  or  equivalently  for  b  =  4. 
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Figure  5.11:  Plot  of  me an  square  error  between  the  eye. Ip  picture  and  its  reconstructed  version 
versus  normalized  number  of  position  and  amplitude  quantization  bits  for  different  values  of 
grid  size. 
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Figure  5.12:  Reconstructed  images  of  the  eye. Ip  picture  via  the  implicit  sampling  strategy  and 
the  iterative  algorithm.  The  number  of  thresholds  is  increased  in  the  sequence  4,6,8  from  top 
to  bottom,  and  the  grid  size  is  increased  in  the  sequence  16,32,64,128  from  left  to  right. 


•  The  number  of  thresholds  which  results  in  smallest  number  of  quantization  bits  is  a 
decreasing  function  of  the  mse.  For  the  curve  corresponding  to  b  =  4,  the  “optimal” 
number  of  thresholds  for  mse  =  1,  .1,  01  are  6,32,  and  128  respectively. 

•  Beyond  a  certain  point,  an  increase  in  the  number  of  thresholds  does  not  result  in  fur¬ 

ther  decrease  in  the  mse,  but  merely  increases  the  total  number  of  bits.  For  the  curves 
corresponding  to  b  =  5,6,  this  phenomenon  is  reached  at  n*  =  16,8  respectively. 

Thus,  the  larger  the  grid  size  is,  the  smaller  the  number  of  thresholds  which  result  in  this 
“saturation”  phenomenon. 

Another  way  to  present  the  quantization  results  of  figure  (5.11)  is  to  plot  the  mean  square 
error  versus  normalized  number  of  position  and  amplitude  bits  as  a  function  of  the  number  of 
thresholds.  This  is  shown  in  the  curves  of  figure  (5.13),  which  correspond  to  reconstruction  from 
a  different  number  of  thresholds.  Various  points  on  each  curve  correspond  to  reconstruction 
with  different  values  of  grid  size.  As  we  would  expect,  the  slopes  of  the  curves  are  negative, 
indicating  that  the  quality  of  reconstruction  improves  as  the  quantization  grid  becomes  finer. 
In  addition,  the  number  of  thresholds  which  results  in  smallest  number  of  quantization  bits  is 
a  function  of  mse.  For  instance,  if  we  are  interested  in  reconstructed  signals  with  mse  <  .556, 
then  the  optimal  number  of  thresholds  is  between  8  and  16.  Finally,  comparing  figures  (5.13) 
and  (5.9),  it  seems  that  for  fixed  quality  of  reconstruction  via  iterative  algorithms,  implicit 
sampling  results  in  lower  number  of  bits  than  semi-implicit  sampling  with  lines  of  unit  slope. 

5.3  Relationship  to  Nyquist  Sampling 

As  we  mentioned  in  Chapter  1,  explicit  sampling  refers  to  schemes  in  which  a  function  is 


represented  by  its  samples  or  derivatives  at  prespecified  points.  An  example  of  such  a  technique 
is  Nyquist  sampling  where  the  amplitude  of  a  signal  is  given  at  equally  spaced  points.  Thus, 
for  a  two-dimensional  BLP  signal  with  N  x  N  region  of  support  in  the  Fourier  domain,  Nyquist 
sampling  involves  amplitude  specification  at  the  nodes  of  a  N  x  N  grid.  Our  main  goal  in  this 
section  is  to  explore  the  relationship  between  Nyquist  sampling/reconstruction  and  reconstruc¬ 
tion  from  level  crossings  via  implicit  or  semi-implicit  sampling  strategies.  We  shall  begin  with 
the  implicit  sampling  strategy  of  Chapter  4. 

5.3.1  Implicit  Sampling 

Recall  from  section  (5.1.1)  that  reconstruction  of  an  N  x  N  signal  from  implicit  samples  of 
its  nt  level  crossings  via  the  linear  least-squares  approach  consists  of  the  following  steps: 

1.  Find  the  level  crossing  contours  associated  with  nt  thresholds. 

2.  Quantize  the  threshold  contours  to  6  bits  by 

•  Superimposing  a  24  x  26  grid  over  the  signal  in  the  space  domain. 

•  Assigning  to  the  center  point  of  all  the  ^  x  ^  squares  which  contain  some  piece 
of  one  or  more  threshold  contours,  the  value  of  the  threshold  whose  contour  comes 
closest  to  the  center  of  the  square. 

3.  Choose  M  of  the  above  ^  x  ^  squares  with  their  center  points  becoming  our  quantized 
reconstruction  samples. 

4.  Reconstruct  the  signal  from  its  quantized  samples  by  finding  its  Fourier  series  coefficients. 

The  above  process  is  shown  pictorially  in  figure  (5.14a).  As  seen  in  figure  (5.14b),  for 
sufficiently  large  number  of  thresholds,  all  the  22i>  squares  associated  with  the  2k  x  26  grid 


will  have  a  threshold  value  associated  with  them.  Moreover,  if  the  grid  size  is  the  same  as 
the  signal  size  i.e.  2*  =  N,  the  number  of  reconstruction  samples,  M,  becomes  equal  to 
N 2,  and  we  will  end  up  with  quantized  samples  of  our  signal  on  a  N  x  N  grid.  Under  these 
circumstances,  the  position  of  our  sampling  points  are  identical  to  that  of  Nyquist  samples,  and 
their  amplitude  corresponds  to  the  threshold  associated  with  them.  Thus,  there  seems  to  be  a 
duality  between  the  above  sampling  set  and  log2  bit  amplitude  quantized  Nyquist  samples. 
However,  it  is  important  to  bear  in  mind  that  these  two  sampling  sets  are  inherently  different 
from  each  other.  In  Nyquist  sampling,  the  amplitude  information  is  quantized  at  the  nodes  of  a 
N  x  N  grid,  whereas,  in  reconstruction  from  threshold  crossings,  the  nodes  of  this  N  x  N  grid 
correspond  to  position  quantized  samples  of  level  crossings.  In  other  words,  Nyquist  sampling 
corresponds  to  amplitude  quantization  at  prespecified  points,  and  reconstruction  from  level 
crossings  corresponds  to  position  quantization  at  prespecified  amplitudes. 

We  will  now  turn  our  attention  to  the  iterative  algorithm  of  section  (4.2.1).  Recall  from 
section  (5.2.2)  that  the  quantization  procedure  for  iterative  reconstruction  of  an  N  x  N  signal 
from  nt  level  crossings  via  the  implicit  sampling  strategy  consists  of  the  following  steps: 

1.  Find  all  the  level  crossing  contours  associated  with  nt  thresholds. 

2.  Quantize  the  threshold  contours  to  b  bits  by  finding  nodes  of  a  2b  x  2b  grid  corresponding 
to  the  boundary  points  of  the  contours. 

3.  Derive  the  range  of  the  signal  for  each  of  the  nodes  of  a  2b  x  2h  grid  by  utilizing  the 
quantized  threshold  contours  of  the  previous  step.  This  will  be  used  as  the  space  domain 
constraint  of  the  algorithm. 

4.  Impose  the  bandlimited  constraint  and  the  space  domain  constraint  of  the  previous  step 


Figure  5.14:  (a)  Position  quantized  samples  obtained  via  im¬ 
plicit  sampling;  (b)  Position  quantized  samples  occupy  all  the 
N2  nodes  of  a  N  x  N  grid  if  the  number  of  thresholds  are 
large  enough. 
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iteratively. 

Thus,  the  iterative  algorithm  derives  the  amplitude  range  for  the  nodes  of  a  2*  x  2*  grid 
by  utilizing  the  quantized  threshold  contours.  Therefore,  the  input  to  the  iterative  algorithm 
can  be  either  represented  by  quantized  threshold  contours,  or  by  the  space  domain  constraint 
obtained  from  them.  Indeed,  as  we  saw  in  section  (5.2.2),  representation  via  the  former  strategy 
involves  tracing  the  boundary  of  threshold  contours,  and  representation  via  the  space  domain 
constraint  involves  specifying  the  range  information  for  the  nodes  of  a  2*  x  2h  grid.  In  the 
case  where  ail  the  contours  associated  with  all  the  thresholds  are  used  for  reconstruction,  the 
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•  points  with  intensity  <  t 

•  points  with  intensity  >  t  and  <  t2 
°  points  with  Intensity  >  1 2 


Figure  5.15:  (a)  Contours  corresponding  to  two  thresholds  tj 
and  t2;  (b)  Amplitude  quantized  samples  derived  from  quan¬ 
tized  threshold  contours  are  identical  to  amplitude  quantized 
Nyquist  samples. 
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intensity  of  each  node  of  the  grid  lies  in  one  of  the  nt  +  1  intervals  defined  by  the  nt  thresholds. 
Thus,  as  shown  in  figure  (5.15),  we  can  think  of  the  nodes  of  the  grid  being  amplitude  quantized 
to  log2(nt  +  1)  bits.  In  addition,  if  the  size  of  the  grid,  26,  assumes  its  minimum  possible  value 

i 

i.e.  N,  then  our  sampling  set  becomes  identical  to  log2(n(  +  1)  bit  amplitude  quantized  Nyquist 
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5.3.2  Semi-implicit  Sampling 

The  discussion  for  the  semi-implicit  sampling  is  somewhat  similar  to  implicit  sampling. 
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Recall  from  section  (5.1.1)  that  reconstruction  of  an  N  x  N  signal  from  its  crossings  with 
nt  arbitrary  functions,  via  the  semi-implicit  sampling  strategy  of  Chapter  2,  consists  of  the 
following  steps: 

1.  Find  all  the  crossings  of  the  signal  with  the  n*  crossing  functions. 

2.  Sample  the  crossings  along  lines  of  rational  slope. 

3.  Choose  M  samples  from  the  intersections  of  the  sampling  lines  and  function  crossings  in 
such  a  way  that  conditions  of  Corollary  (2.3)  are  satisfied. 

4.  Quantize  the  position  of  the  chosen  samples  to  b  bits. 

5.  Reconstruct  the  signal  from  its  quantized  samples  by  finding  its  Fourier  series  coefficients. 

The  above  process  is  shown  pictorially  in  figure  (5.16a).  Suppose  that  the  following  conditions 
are  satisfied: 

•  The  crossing  functions  are  chosen  to  be  constants,  so  that  function  crossings  become  level 
crossings. 

•  The  sampling  lines  are  chosen  to  be  horizontal  or  vertical  and  equally  spaced. 

•  The  number  of  reconstruction  samples  is  IV2. 

•  The  number  of  quantization  bits  is  log2  N . 

As  seen  in  figure  (5.16b),  under  these  circumstances,  the  location  of  our  N2  position  quantized 
samples  corresponds  to  the  nodes  of  a  N  x  N  grid,  or  equivalently  the  location  of  Nyquist 
samples,  and  their  amplitudes  correspond  to  the  threshold  crossing  associated  with  them.  Note 
that  for  the  above  conditions  to  hold,  the  number  of  thresholds  used  for  reconstruction  must 
be  large  enough  so  that  each  of  the  N  horizontal  (or  vertical)  lines  contain  N  samples  after 

149 


\ 

i 

» 


» 


WlWUWJUUWWIUljUIIOTWW  W  W*V*  1 JNJWTC  HR  «  PJI W  JUl  PJI  Ml  lU^WWWWWGWW 


quantizing  the  i  (or  y)  position  of  the  samples  to  log2  N  bits.  Thus,  similar  to  the  implicit 
sampling  case,  there  seems  to  be  a  duality  between  the  samples  obtained  in  the  above  manner 
and  log2  nt  bit  amplitude  quantized  Nyquist  samples.  However,  it  is  important  to  bear  in  mind 
that  these  sampling  sets  are  inherently  different  from  each  other.  In  Nyquist  sampling,  the 
amplitude  information  is  quantized  at  the  nodes  of  a  JV  x  JV  grid,  whereas,  in  semi-implicit 
reconstruction  from  threshold  crossings,  the  nodes  of  this  N  x  N  grid  correspond  to  samples  of 
level  crossings  position  quantized  along  the  sampling  lines.  In  other  words,  Nyquist  sampling 
corresponds  to  amplitude  quantization  at  prespecified  points,  and  reconstruction  from  level 
crossings  corresponds  to  position  quantization  along  sampling  lines  at  prespecified  amplitudes. 

We  will  now  turn  our  attention  to  the  iterative  algorithm  of  section  (3.3.1).  As  shown  in 
figure  (5.17a)  ,  the  iterative  algorithm  utilizes  the  position  quantized  semi-implicit  samples  in 
order  to  derive  the  space  domain  constraint.  Therefore,  the  input  to  the  iterative  algorithm 
can  be  represented  either  by  position  quantized  samples,  i.e.  in  a  similar  fashion  to  the  linear 
least-squares  and  recursive  approach,  or  by  the  intensity  range  of  equally  spaced  points  on  the 
sampling  lines.  In  the  latter  case,  since  all  the  intersections  of  sampling  lines  with  the  threshold 
contours  are  utilized,  the  intensity  of  the  equally  spaced  points  lie  in  one  of  nt  +  1  intervals 
defined  by  the  r»t  thresholds.  In  addition,  suppose  that  the  following  conditions  are  satisfied: 

•  We  have  JV  equally  spaced  horizontal  or  vertical  sampling  lines. 

•  The  number  of  equally  spaced  points  on  sampling  lines,  M,  is  equal  to  JV. 

As  figure  (5.17b)  shows,  under  these  circumstances,  the  JV  equally  spaced  samples  on  JV  hori¬ 
zontal  or  vertical  lines  correspond  to  nodes  of  a  JV  x  JV  grid,  and  since  the  amplitude  of  each 


•  Non-quantlzed  samples 
o  Quantized  samples 


Sampling  lines 


Lines  used  for  quantization 
along  sampling  lines 


(a)  (b) 

Figure  5.16:  (a)  Position  quantized  samples  obtained  via 
semi-implicit  sampling;  (b)  If  we  have  N  equidistant  horizon¬ 
tal  or  vertical  sampling  lines,  and  the  number  of  thresholds 
is  large  enough,  the  position  quantized  samples  will  occupy 
all  the  N2  nodes  of  a  N  x  N  grid. 


identical  to  log2(nt  +  1)  bit  amplitude  quantized  Nyquist  samples. 

To  summarize  this  section,  we  have  found  conditions  under  which  the  location  of  position 
quantized  samples  obtained  via  semi-implicit  and  implicit  sampling  become  identical  to  those 
of  Nyquist  samples.  In  addition,  under  certain  circumstances,  the  range  information  used  as 
an  input  to  the  iterative  algorithms  can  be  considered  to  be  identical  to  amplitude  quantized 
Nyquist  samples.  Thus,  it  appears  that  representation  of  two-dimensional  signals  via  their 
amplitude  quantized  explicit  samples  is  intimately  related  to  their  position  quantized  semi- 
implicit  or  implicit  samples,  and  that  reconstruction  from  multiple  level  threshold  crossings  has 
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•  points  with  Intensity  <  t, 
o  points  with  intensity  >  1 1  and  <  t2 

o  points  with  Intensity  >  t2 


Figure  5.17:  (a)  Derivation  of  space  domain  constraint  for  the 
iterative  algorithm  of  semi-implicit  sampling;  (b)  If  we  have 
N  equidistant  horizontal  or  vertical  sampling  lines,  together 
with  N  equally  spaced  samples  on  each  line,  our  sampling  set 
corresponds  to  amplitude  quantized  Nyquist  samples. 


bridged  the  gap  between  explicit,  semi-implicit  and  implicit  sampling  strategies. 


5.4  Summary 


In  this  chapter,  we  have  presented  a  preliminary  investigation  of  quantization  properties  of 
various  sampling  and  reconstruction  schemes  as  a  function  of  the  number  of  thresholds.  We 
started  with  the  linear  least-squares  approach  via  the  QR  decomposition,  which  is  our  most 
stable  reconstruction  algorithm  for  both  the  semi-implicit  and  implicit  sampling  strategies. 


There  are  a  variety  of  factors  influencing  the  robustness  of  reconstruction  via  the  linear  least- 
squares  approach.  In  section  (5.1.2),  we  found  that  for  fixed  number  of  thresholds  and  a 
given  sampling  strategy,  the  number  of  position  bits  is  initially  decreased  as  the  number  of 
reconstruction  samples  is  increased  from  its  minimum  value.  However,  as  seen  in  figure  (5.4), 
increasing  the  number  of  samples  beyond  a  certain  point  will  merely  result  in  an  increase 
in  total  number  of  bits  used  for  representing  the  signal.  In  section  (5.1.3),  we  found  that 
for  fixed  number  of  reconstruction  samples,  the  quantization  characteristics  initially  improve 
as  the  number  of  thresholds  is  increased.  However,  for  large  number  of  thresholds,  further 
increase  in  the  number  of  thresholds  does  not  necessarily  lead  to  more  favorable  quantization 
characteristics. 

As  far  as  iterative  algorithms  go,  the  “optimum”  number  of  thresholds,  which  results  in 
smallest  number  of  position  and  amplitude  bits,  is  highly  dependent  on  the  quality  of  the 
reconstructed  images.  More  specifically,  as  seen  in  figures  (5.9)  and  (5.11),  for  smaller  values 
of  the  mean  square  error,  the  “optimum”  number  of  thresholds  is  larger. 

Recall  that  our  goal  in  this  chapter  has  been  to  demonstrate  that  the  quantization  charac¬ 
teristics  of  our  sampling  and  reconstruction  algorithms  lie  in  between  Nyquist  and  zero  crossing 
sampling  and  reconstruction.  While  recovery  of  an  IV  x  IV  signal  from  its  Nyquist  samples  re¬ 
quires  minimum  number  of  position  bits,  i.e.  log2  N ,  and  maximum  number  of  amplitude  bits, 
i.e.  infinite,  and  recovery  from  zero  crossings  requires  maximum  number  of  position  bits,  i.e. 
infinite,  and  minimum  number  of  amplitude  bits,  i.e.  1,  the  position  and  amplitude  quantiza¬ 
tion  requirements  of  our  sampling  strategies  for  reconstruction  from  multiple  level  crossings  lie 
in  between  these  two  extremes.  In  fact,  the  experimental  results  of  this  chapter  seem  to  suggest 
that  the  optimal  number  of  thresholds,  which  results  in  minimum  number  of  total  amplitude 


tion  requirements  of  our  sampling  strategies  for  reconstruction  from  multiple  level  crossings  lie 
in  between  these  two  extremes.  In  fact,  the  experimental  results  of  this  chapter  seem  to  suggest 
that  the  optimal  number  of  thresholds,  which  results  in  minimum  number  of  total  amplitude 
and  position  bits,  is  neither  infinite,  as  it  is  with  Nyquist  sampling,  nor  is  it  one,  as  is  the  case 
with  zero  crossing  sampling.  Indeed,  this  optimum  number  depends  on  a  variety  of  factors  such 


Quality  of  reconstruction. 


The  specific  sampling  strategy  used  i.e.  semi-implicit  or  implicit  sampling. 


The  specific  reconstruction  strategy. 


The  number  of  reconstruction  samples. 


Finally,  as  we  saw  in  section  (5.3),  representation  of  two-dimensional  signals  via  their  am¬ 
plitude  quantized,  explicit,  Nyquist  samples  is  intimately  related  to  their  position  quantized 
implicit  or  semi-implicit  samples.  The  preliminary  results  of  this  chapter  seem  to  indicate  that 
reconstruction  from  multiple  level  threshold  crossings,  has  bridged  the  gap  between  explicit, 
semi-implicit,  and  implicit  sampling  strategies,  unified  seemingly  unrelated  sampling  schemes, 
and  provided  us  with  a  spectrum  of  sampling  techniques  for  multidimensional  signals. 
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Chapter  6 


Summary,  Conclusions,  and  Future 
Directions  of  Research 


Our  main  goal  in  this  thesis  has  been  to  derive  sampling  and  reconstruction  schemes  whose 
characteristics  lie  in  between  Nyquist  sampling  and  zero  crossing  sampling  as  proposed  by  Curtis 
[51].  As  we  mentioned  in  Chapter  1,  while  the  required  bandwidth  for  representing  signals  via 
their  Nyquist  samples  is  minimal,  to  be  able  to  recover  the  signal  exactly,  the  dynamic  range 
or  equivalently  the  number  of  amplitude  bits  must  be  infinite.  On  the  other  hand,  since  in 
representing  multidimensional  signals  via  their  one  level  threshold  crossings  the  location  of  the 
level  crossings  must  be  known  with  large  accuracy,  the  required  bandwidth  is  extremely  large, 
while  the  dynamic  range  is  minimum.  Thus  our  objective  has  been  to  derive  sampling  strategies 
whose  bandwidth  and  dynamic  range  characteristics  lie  in  between  these  two  extremes,  and  can 
be  used  for  reconstruction  of  signals  from  their  multiple  level  crossings. 

Our  approach  to  solving  this  problem  has  been  to  formulate  it  in  terms  of  multivariate 
interpolation  theory.  The  reasons  behind  this  formulation  are  two-fold.  First,  BLP  signals  which 
encompass  a  fairly  large,  general  class  of  signals,  can  be  represented  in  terms  of  polynomials. 
Second,  a  rather  large  body  of  the  mathematical  results  in  polynomial  interpolation  theory  can 
be  applied  directly  to  the  problem  at  hand.  As  we  mentioned  in  Chapter  2,  although  univariate 


polynomial  interpolation  is  relatively  straightforward,  the  inherent  difficulty  in  multivariate 
interpolation  is  the  fact  that  unlike  the  univariate  case,  there  are  no  Chebychef  systems  in  R ‘ 
for  s  >  2.  Indeed,  our  main  theoretical  results  in  this  thesis  deal  with  two  major  strategies  for 
overcoming  this  difficulty. 

Our  first  strategy,  described  in  Chapter  2,  consisted  of  imposing  restrictions  on  the  location 
of  interpolation  points,  used  for  recovery  of  the  bivariate  polynomial  associated  with  the  signal 
under  consideration.  To  this  end,  we  developed  several  theoretical  results  on  multivariate 
polynomial  interpolation  theory  using  algebraic  geometric  concepts.  We  then  used  these  results 
to  derive  a  variety  of  semi-implicit  sampling  strategies  to  provide  us  with  sufficient  conditions 
under  which  multidimensional  BLP  signals  can  be  recovered  from  their  non-uniform  samples 
on  lines  of  rational  slope.  To  utilize  these  results  in  the  context  of  reconstruction  from  level 
crossings,  the  non-uniform  samples  were  chosen  at  the  intersection  of  sampling  lines  with  level 
crossings.  As  we  saw  in  section  (2.3),  the  semi-implicit  sampling  strategy  can  also  be  applied 
to  a  variety  of  other  problems  such  as  reconstruction  from  crossings  of  the  signal  with  arbitrary 
functions,  and  reconstruction  of  multidimensional  signals  from  their  projections. 

The  major  drawback  of  the  line  sampling  strategy  for  reconstruction  from  level  crossings 
is  the  fact  that  in  general,  we  are  never  guaranteed  to  get  enough  intersections  between  the 
sampling  lines  and  level  crossings  to  satisfy  the  conditions  of  the  theorems  in  Chapter  2.  As 
we  saw  in  section  (3.4),  although  we  can  improve  the  likelihood  of  satisfying  these  conditions 
by  careful  choice  of  the  slope  and  position  of  sampling  lines,  the  basic  problem  still  remains  in 
the  sense  that  the  semi-implicit  sampling  strategy  is  incapable  of  reconstructing  signals  from 
an  arbitrary  number  of  thresholds. 

To  overcome  the  above  difficulty,  in  Chapter  4,  we  proposed  a  second  approach  to  the 
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problem  of  reconstruction  from  level  crossings.  While  the  semi-implicit  sampling  strategy  of 
Chapter  2,  can  be  applied  to  the  more  general  problem  of  reconstruction  from  non-uniformly 
spaced  samples,  the  implicit  sampling  strategy  of  Chapter  4  can  only  be  applied  to  reconstruc¬ 
tion  from  level  crossings  or  crossings  with  functions,  whose  bandwidth  lies  within  the  bandwidth 
of  the  signal.  The  main  advantage  of  the  implicit  sampling  scheme,  however,  is  that  unlike  the 
semi-implicit  approach  it  can  reconstruct  signals  from  an  arbitrary  number  of  thresholds.  The 
major  result  in  Chapter  4  states  that  for  almost  all  signals  with  N  x  N  region  of  support  in  the 
Fourier  domain,  almost  all  A:  >  0  points  from  its  a  level  crossings  and  IV2  -  k  points  from  its  /? 
level  crossings  are  sufficient  to  uniquely  specify  it.  This  result  was  extended  to  situations  where 
the  number  of  thresholds  is  larger  than  2  and  to  the  problem  of  reconstruction  from  crossings 
with  sinusoids,  whose  frequencies  lie  within  the  bandwidth  of  the  signal. 

Having  developed  the  semi-implicit  and  implicit  sampling  strategies,  we  then  proposed  a 
variety  of  reconstruction  algorithms  in  Chapters  3  and  4.  The  most  straightforward  way  of 
carrying  out  reconstructions  for  both  the  semi-implicit  and  implicit  sampling  strategies,  is  to 
solve  a  linear  system  of  equations  to  find  the  Fourier  series  coefficients  associated  with  the  signal. 
We  proposed  two  algorithms  for  solving  the  linear  least-squares  problem:  QR  decomposition 
and  the  conjugate  gradient  algorithm.  While  QR  decomposition  is  considerably  more  robust 
than  the  conjugate  gradient  algorithm,  its  storage  requirements  for  a  signal  with  N  x  N  region  of 
support  in  the  Fourier  domain  are  of  the  order  of  N4.  We  also  proposed  the  recursive  algorithm 
of  section  (3.2)  in  conjunction  with  the  semi-implicit  sampling  strategy;  the  recursive  approach 
is  not  storage  intensive  and  is  efficient  computationally,  although  it  is  somewhat  ill  conditioned. 
To  circumvent  the  problems  associated  with  the  recursive  and  linear  least-squares  approach, 
we  proposed  the  line  by  line  iterative  algorithm  of  section  (3.3.1)  for  the  semi-implicit  sampling 


strategy,  and  the  iterative  algorithm  of  section  (4.2)  for  the  implicit  sampling  strategy.  As  it 
turns  out,  these  iterative  algorithms  are  relatively  stable  and  their  storage  and  computational 
requirements  are  not  demanding. 

The  input  requirements  of  the  iterative  algorithms  however,  are  somewhat  different  from 
our  other  algorithms,  since  they  utilize  all  the  available  level  crossing  information  to  derive  the 
space  domain  constraint  for  the  iterations.  More  specifically,  the  iterative  algorithm  of  section 
(3.3.1)  requires  all  the  intersections  of  sampling  lines  and  threshold  contours,  and  the  iterative 
algorithm  of  section  (4.2)  requires  all  the  quantized  threshold  contours.  An  important  feature 
of  the  iterative  algorithm  of  section  (3.3.1)  is  the  fact  that,  if  the  number  of  threshold  crossings 
on  each  sampling  line  exceeds  the  number  of  Fourier  harmonics  of  the  one-dimensional  signal 
associated  with  the  line,  then  the  solution  obtained  via  the  iterative  algorithm  is  unique,  only 
in  the  limit  as  the  number  of  equally  spaced  samples  on  the  sampling  lines,  A/,  tends  to  infinity. 
Since  in  practice  M  can  only  be  finite,  the  solution  obtained  via  the  iterative  algorithm  is  only 
an  approximate  one.  As  far  as  the  implicit  sampling  strategy  goes,  while  the  linear  least-squares 
approach  requires  N 2  samples  of  threshold  crossings  in  order  to  recover  an  N  x  N  signed,  the 
iterative  algorithm  requires  quantized  versions  of  all  the  threshold  contours  themselves.  Similar 
to  the  semi-implicit  case,  the  solution  obtained  via  the  iterative  algorithm  is  unique  only  in 
the  limit  as  the  size  of  quantization  grid  approaches  zero.  Thus  the  solution  obtained  via  the 
iterative  algorithm  with  finite  grid  size  is  only  an  approximate  one. 

In  Chapter  5,  we  presented  a  preliminary  investigation  of  the  quantization  properties  of 
some  of  our  reconstruction  algorithms  as  a  function  of  the  number  of  thresholds.  We  started 
with  the  linear  least-squares  approach  via  the  QR  decomposition,  which  is  our  most  stable 
reconstruction  algorithm  for  both  the  semi-implicit  and  implicit  sampling  strategies.  As  we 
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saw,  there  are  a  variety  of  factors  influencing  the  robustness  of  reconstruction  via  the  linear 
least-squares  approach.  In  section  (5.1.2),  we  found  that  for  a  fixed  number  of  thresholds  and 
a  given  sampling  strategy,  the  number  of  position  bits  is  initially  decreased  as  the  number  of 
reconstruction  samples  is  increased  from  its  minimum  value.  However,  as  seen  in  figure  (5.4), 
increasing  the  number  of  samples  beyond  a  certain  point  will  merely  result  in  an  increase  in 
total  number  of  bits  used  for  representing  the  signal.  In  section  (5.1.3),  we  found  that  for  a  fixed 
number  of  reconstruction  samples,  the  quantization  characteristics  of  the  linear  least-squares 
approach  initially  improve  as  the  number  of  thresholds  is  increased.  However,  for  a  large  number 
of  thresholds,  further  increase  in  the  number  of  thresholds  does  not  necessarily  lead  to  more 
favorable  quantization  characteristics.  As  we  saw  in  section  (5.2),  the  quantization  properties 
of  the  iterative  algorithms  suggest  that  the  “optimum”  number  of  thresholds  which  results  in 
the  smallest  number  of  position  and  amplitude  bits  is  highly  dependent  on  the  quality  of  the 
reconstructed  images.  More  specifically,  as  figures  (5.9)  and  (5.11)  show,  for  smaller  values 
of  the  mean  square  error,  the  “optimum”  number  of  thresholds  is  larger.  It  is  important  to 
mention  that  the  results  presented  in  Chapter  5  are  rather  preliminary  and  that  the  conclusions 
are  somewhat  tentative.  Our  hope  is  that  these  speculative  results  can  be  used  as  a  starting 
point  for  further  research  in  the  applications  of  the  theory  to  areas  of  multidimensional  signal 
representation  and  image  coding. 

Recall  that  our  goal  in  this  thesis  has  been  to  derive  sampling  schemes  whose  bandwidth 
and  dynamic  range  characteristics  lie  in  between  those  of  Nyquist  and  zero  crossings  sampling. 
While  recovery  of  an  N  x  N  signal  from  its  Nyquist  samples  requires  minimum  number  of 
position  bits,  i.e.  log2N,  and  maximum  number  of  amplitude  bits,  i.e.  infinite,  and  recovery 
from  zero  crossings  requires  maximum  number  of  position  bits,  i.e.  infinite  and  minimum 


number  of  amplitude  bits,  i.e.  1,  the  position  and  amplitude  quantization  requirements  of  our 
sampling  strategies  for  reconstruction  from  multiple  level  crossings  lie  in  between  these  two 
extremes.  In  fact,  the  experimental  results  of  Chapter  5  seem  to  suggest  that  the  optimal 
number  of  thresholds,  which  results  in  minimum  number  of  total  amplitude  and  position  bits, 
is  neither  infinite,  as  it  is  with  Nyquist  sampling,  nor  is  it  one,  as  is  the  case  with  zero  crossing 
sampling.  Indeed,  this  optimum  number  depends  on  a  variety  of  factors  such  as: 

•  Quality  of  reconstruction. 

•  The  specific  sampling  strategy  used  i.e.  semi-implicit  or  implicit  sampling. 

•  The  specific  reconstruction  strategy. 

•  The  number  of  reconstruction  samples. 

Finally,  as  we  saw  in  section  (5.3),  representation  of  two-dimensional  signals  via  their  am¬ 
plitude  quantized  explicit  samples  is  intimately  related  to  their  position  quantized  implicit  or 
semi-implicit  samples.  The  results  of  Chapter  5  seem  to  indicate  that  not  only  does  the  am¬ 
plitude  and  position  quantization  characteristics  of  our  sampling  and  reconstruction  schemes 
lie  in  between  those  of  Nyquist  and  zero-crossings,  but  also  under  certain  circumstances,  semi- 
implicit  and  implicit  sampling  strategies  become  a  special  case  of  Nyquist  sampling.  In  short, 
reconstruction  from  multiple  level  threshold  crossings,  has  bridged  the  gap  between  explicit, 
semi-implicit  and  implicit  sampling  strategies,  unified  seemingly  unrelated  sampling  schemes, 
and  provided  us  with  a  spectrum  of  sampling  techniques  for  multidimensional  signals. 


6.1  Directions  for  Future  Research 


A  wide  variety  of  different  directions  are  possible  for  expanding  or  extending  the  results 
of  this  thesis.  The  most  natural  application  of  the  results  are  in  the  area  of  multidimensional 
signal  representation  and  coding.  As  we  mentioned  earlier,  the  quantization  results  of  Chapter 
5  are  rather  preliminary,  and  most  of  the  conclusions  are  somewhat  tentative.  However,  these 
speculative  results  can  be  used  as  a  starting  point  for  a  rigorous  and  thorough  investigation  of 
the  quantization  characteristics  of  the  various  sampling  and  reconstruction  schemes  presented  in 
this  thesis.  Such  characterization  involves  consideration  of  many  coding  issues,  a  large  number 
of  experiments  and  comparison  to  other  existing  coding  schemes. 

Our  theoretical  results  can  also  be  applied  to  the  area  of  multidimensional  signal  reconstruc¬ 
tion  from  projections.  This  problem  arises  in  diverse  fields  as  X-ray  tomography,  transmission 
electron  microscopy,  and  radio  astronomy.  In  section  (2.3.2),  we  used  Theorem  (2.7)  to  gener¬ 
alize  the  one-projection  theorem  due  to  Mersereau  and  Oppenheim  [45]  to  a  multiple-projection 
theorem.  The  major  problem  in  reconstructing  bandlimited  functions  of  order  M  from  a  single 
projection  is  a  computational  one,  and  is  due  to  the  high  order  of  the  polynomials  involved. 
More  specifically,  in  order  to  recover  all  M2  samples  from  a  single  projection,  it  is  necessary  to 
work  with  a  polynomial  of  degree  greater  than  or  equal  to  M2 .  However,  reconstruction  from 
multiple  projections  at  angles  specified  by  Corollary  (2.5),  results  in  lower  degree  polynomials 
and  thus,  a  more  stable  numerical  problem.  Extensive  simulations  are  needed  to  verify  this 
idea  experimentally,  and  to  compare  it  with  the  existing  reconstruction  schemes. 

Our  results  in  bivariate  polynomial  interpolation  can  also  be  applied  to  a  variety  of  other 
problems.  An  example  of  such  problems  is  multidimensional  FIR  filter  design,  where  the  fre- 


quency  response  of  the  filter  is  expressed  in  terms  of  a  two-dimensional  polynomial.  If  the 
frequency  response  of  the  filter  is  specified  at  points  located  on  lines  of  rational  slope,  we  can 
apply  a  variety  of  the  techniques  described  in  Chapter  3  to  determine  the  impulse  response  of 
the  filter.  Another  application  of  the  multivariate  polynomial  interpolation  is  surface  interpo¬ 
lation  from  scattered  data.  This  problem  arises  frequently  in  fields  such  as  geology,  astronomy, 
oceanography,  and  machine  vision  where  it  is  necessary  to  recover  a  signal  or  function  whose 
values  are  known  along  certain  contours,  paths,  or  curves.  For  instance  in  machine  vision,  pri¬ 
mal  sketch  descriptions  of  several  images  are  matched,  either  by  stero  or  motion  computation, 
to  obtain  a  description  of  surface  information  at  the  zero  crossings  contours  of  the  convolution 
of  the  image  with  the  Laplacian  of  a  Gaussian. 

As  we  saw  in  section  (2.3.1),  the  theoretical  results  of  Chapter  2  can  be  used  to  reconstruct 
signals  from  their  crossings  with  arbitrary  functions.  A  potential  application  of  this  result  is  in 
the  conversion  of  halftone  images  to  contone  ones.  The  halftone  process  has  been  used  for  more 
than  a  century  to  convert  continuous  tone  pictures  into  a  regular  patterns  of  black  and  white 
dots  which  can  then  be  printed.  Mathematically  speaking,  the  halftone  version  of  a  continuous 
tone  image  can  be  obtained  by  comparing  the  value  of  the  signal  with  a  two-dimensional  periodic 
function,  called  the  screen  function,  and  producing  a  white  or  black  pixel  on  a  high  contrast 
medium  depending  on  whether  the  signal  value  is  higher  or  lower  than  the  screen  function.  In 
effect,  the  boundary  of  black  and  white  pixels  in  the  halftone  image  corresponds  to  the  crossing 
of  the  signal  with  the  screen  function.  Thus,  the  theoretical  results  of  Chapter  2  can  be  applied 
to  reconstruct  the  continuous  signal  from  its  crossings.  Few  examples  of  such  reconstructions  via 
the  conjugate  gradient  and  recursive  algorithms  were  shown  in  Chapter  3.  In  these  examples, 
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the  screen  function  was  chosen  to  be 

A  {1  +  cos[2jr(pi  +  gy)]} 

Thus,  future  work  must  focus  on  the  applicability  of  these  results  to  more  general  screen  func¬ 
tions,  numerical  robustness  of  reconstruction  in  practical  situations,  and  comparison  with  ex¬ 
isting  conversion  techniques. 

As  we  mentioned  in  Chapter  1,  a  major  drawback  of  reconstruction  from  one-level  crossings 
is  that  the  location  of  the  crossings  need  to  be  known  extremely  accurately.  This  accuracy  lim¬ 
its  the  applicability  of  reconstruction  from  one-level  crossings  in  many  practical  situations  such 
as  image  restoration.  More  specifically,  consider  an  image  which  has  undergone  nonintentional 
nonlinear  processing  such  as  a  high  contrast  recording  medium  with  few  intensity  levels.  If  it  is 
distorted  in  such  a  way  as  to  preserve  one  or  more  of  its  level  crossings,  it  is  possible,  at  least 
in  principle,  to  recover  it  from  the  level  crossings  of  its  distorted  version.  Since  reconstruction 
from  a  single  threshold  needs  very  accurate  location  of  the  one-level  crossings,  successful  recov¬ 
ery  of  the  original  image  requires  an  extremely  high  resolution  recording  medium.  However, 
considering  the  results  of  Chapter  5,  we  would  expect  this  requirement  to  be  substantially  re¬ 
duced  if  the  original  image  is  reconstructed  from  multiple  level  crossings  of  its  distorted  version. 
Thus,  the  spectrum  of  sampling/reconstruction  techniques  described  in  this  thesis  enables  us  to 
store/retrieve  data  from  a  wide  variety  of  recording  media  with  different  resolution  and  dynamic 
range  characteristics.  Similar  remarks  can  be  made  about  analog  communication  channels  with 
different  bandwidth  and  dynamic  range  characteristics. 

From  a  theoretical  point  of  view,  we  can  extend  our  results  in  a  number  of  different  ways. 
We  have  not  addressed  the  robustness  of  our  sampling  and  reconstruction  schemes  with  re- 


spect  to  the  bandlimitedness  and  periodicity  assumption.  We  would  expect  the  level  crossing 
contours  of  an  almost  BLP  signal  to  be  somewhat  similar  to  the  threshold  contours  of  its  ban- 
dlimited,  periodic  version.  Further  investigation  is  needed  to  determine  how  violations  of  these 
assumptions  affect  the  quality  of  reconstruction  as  a  function  of  the  number  of  thresholds.  A 
different  problem  involves  finding  sufficient  conditions  for  unique  reconstruction  of  nonperiodic 
signals.  These  conditions  have  already  been  found  for  the  case  of  reconstruction  from  one-level 
crossings  [51], 

Finally,  the  results  in  this  thesis  could  potentially  have  a  major  impact  in  the  area  of  signal 
representation  and  manipulation.  In  most  signal  processing  applications,  multidimensional 
signals  such  as  images  are  usually  represented  by  a  two-dimensional  array  consisting  of  the 
intensity  of  the  signal  on  a  square  or  rectangular  array,  with  most  of  the  processing  done  in  the 
space  or  frequency  domain.  In  this  thesis,  we  have  shown  that  multidimensional  signals  can  be 
robustly  represented  with  an  arbitrary  number  of  threshold  contours  or  their  samples.  Thus, 
one  might  consider  ways  of  processing  signals  in  the  threshold  domain  as  opposed  to  the  more 
conventional  space  or  frequency  domains.  More  specifically,  signal  processing  operations  which 
preserve  the  bandlimitedness  and  periodicity  of  the  signal,  can  be  carried  out  in  the  threshold 
domain  by  simply  operating  on  few  threshold  contours  of  the  signal  in  order  to  obtain  a  new 
set  of  threshold  contours  of  the  “processed”  signal. 

An  interesting  application  of  the  above  idea  might  be  in  the  area  of  morphological  analysis 
of  images  where  some  of  the  morphological  operations  such  as  erosion  and  dialation  commute 
with  thresholding.  One  way  to  carry  out  these  operations  on  a  function  /,  at  least  in  principle, 
is  to  process  all  the  level  crossings  of  /  separately  in  order  to  “build”  the  desired  signal  by 
stacking  its  level  crossings.  Thus,  in  situations  where  we  know  that  the  dialated  or  eroded 
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version  of  a  signal  /  is  bandlimited  and  periodic,  we  can  perform  these  operations  on  few  level 
crossings  of  /  instead  of  all  of  them.  Clearly,  most  of  the  above  ideas  are  rather  speculative, 
and  further  research  is  needed  to  verify  them  both  theoretically  and  experimentally. 
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Appendix  A 

Proof  of  Theorem  3.8 


In  this  appendix  we  will  prove  Theorem  (2.4)  of  section  (2.2.2). 


Theorem  A.l  Let  /o,...,/p  be  distinct  curves  with  U,  the  ith  curve  given  by: 

z  =  a,wm  oii  7*  0 

where  m  <  n  is  an  arbitrary  integer,  p  is  the  smallest  integer  such  that 

P 

^[(m  +  l)n  -  2mi  +  1]  >  (n  +  l)2 

i= o 

and 

{(ti>j**,r^)|j  =  0,  ...,(m+  l)n  -  2mt;}  (A.l) 

the  set  of  distinct  points  on  li.  Jf  none  of  the  interpolation  points  defined  by  (A.l)  are  equal  to 
(0,0),  the  common  intersection  of  all  the  curves,  then  for  any  data  set 

{t^\j  =  0, ...,  (m  +  l)n  -  2m»;  i  =  0,  ...,p;  } 

there  is  a  unique  bivariate  polynomial  of  the  form 

p{w,z)  =  f^a(*’,.7')u;V 
i=o  j=o 

such  that 

p(wj,\«j’))  =  tV  j  =  0,...,(m+  l)n-  2mt;  »' =  0,  ...,n; 


Proof: 

Our  approach  will  be  to  describe  a  method  of  reconstruction  of  the  polynomial  p(w,  z)  which 
satisfies  the  conditions  described  in  the  theorem.  In  the  process  of  reconstruction,  we  will  show 


that  this  polynomial  exists  and  is  unique.  The  reconstruction  algorithm  consists  of  p+  1  steps. 
We  will  use  induction  to  show  that  in  the  tth  step  we  can  find  the  2mi  coefficients  given  by  the 
set 

{a(/i, /2)l^i  +  ^2  =  (>  ~  1)”*,  —,im  -  1,  n(m  +  1)  -  t'm  +  1,  +1)} 

Before  we  start  the  induction,  it  is  worthwhile  to  mention  that  sampling  the  bivariate 
polynomial  p(w,z)  along  the  tth  curve,  ft,  which  is  given  by 

Z  =  OiWm 

is  equivalent  to  sampling  the  univariate  polynomial 

n(m+I) 

p,(w)  =  p(w,a,wm)  =  ^2  br)wT  (A-2) 

r=0 

where 

bll)  =  £  £  ofraUM  (A.3) 

I  j=t»  fo=l» 

I I  +  mi  2 =r 

In  the  first  step  of  the  algorithm,  we  can  use  the  points  on  the  curve  l0  given  by  the  set 

{(wj0),aowj0))|i  =  0,...,(m+  l)n} 
together  with  their  corresponding  data  set 

{tj0,|j  =  °,...,(m+  l)n} 

in  order  to  find  for  i  =  0, ...,  n(m+  1).  This  can  be  done  because  po(w)  of  equation  (A.l)  is 
a  one  dimensional  polynomial  of  degree  n(m  +  1)  and  thus  any  (m  +  l)n  +  1  distinct  samples 
of  it  are  sufficient  to  find  its  coefficients.  We  can  now  use  the  value  of  the  quantities 

{b\0)\i  =  0,  ...,m  -  1,  n(m  4-  1)  -  (m  -  l),...,n(m+  1);  } 

together  with  equation  (A.3)  to  find  the  coefficients 

{<*(/!, /2) |/i  +  ml2  =  0,  ...,m  -  1, n(m  +  1)  -  (m  -  l),...,n(m+  1);  } 

As  we  will  see  later,  the  values  of  the  remaining  coefficients  of  the  polynomial  po(«>),  which 
are  found  in  the  first  step  of  the  algorithm,  will  be  used  in  future  steps.  More  specifically,  for 
j  =  0,  ...,m  -  1  the  quantities  and  will  be  used  in  the  (»  +  l)st  step  of  the 

algorithm. 

Having  shown  the  validity  of  the  induction  hypothesis  for  the  first  step  of  the  algorithm,  we 
will  now  show  that  if  in  steps  1  through  i  the  quantities 

{a(li,li)\li  +  mi2  =  0,  -  1,  n(m  +  1)  -  »m  +  1,  ...,n(m  +  1);  } 

and 

=  0, r  =  jm,jm+  l,...,n(m+  1)  -  jm  +  1;} 


are  found,  then  in  the  (t  +  l)st  step  the  quantities 

-  +mf2  =  mi,  ...,m(i  +  1)  -  1,  n(m  +  1)  -  (t  +  l)m, n(m+  1)  -  im;  } 

and 

{6^,+1|r  =  im, n(m  +  1)  -  t'm;  } 
can  be  evaluated.  Rearranging  the  terms  in  equation  (A.l)  we  get: 

p{w,ciiWm)  =  p(w,OLiWm)~  ^2  X]  a^a(lul2)wl‘+h 

ij=li  /  2=n 

ii+mi2==0,...,m-*l,  n(m+  l)-tm+l,...,n(m+l) 

£  £  c£a{h,l*Wl+la 

i1=ti  i2=  ii 

=  Y,  br)v>r  (A-4) 

r—\m 

By  hypothesis,  since  the  point  (0,0)  is  on  the  intersection  of  all  the  curves  l0,...,lp,  it  could  not 
possibly  be  one  of  the  interpolation  points.  Therefore  we  have 

u/j’*  ^  0  j  =  0,  +  1)  -  2im; 

This  implies  that  the  points  on  the  ith  curve,  given  by  the  set 

=  0,  +  1)  -  2tm} 

are  sufficient  to  uniquely  specify  the  coefficients 

{6^|r  =  im,  ...,n(m  +  1)  -  im;  } 

of  the  univariate  polynomial  given  by  equation  (A. 4).  The  values  of  b's  found  in  the  (i  +  l)st 
step  of  the  algorithm  together  with  the  ones  found  in  previous  steps  can  now  be  used  for  finding 
the  coefficients  of  p(w,z).  More  specifically,  for  j  =  0,  ...,m  —  1  the  quantities 

{bim+j\k  =  o,  •,«';} 

can  be  used  to  find  the  coefficients 

{a{h,h)\h  +  mh  =  mi  +  j,} 

Using  equation  (A. 3)  we  have 

bim+j  =  £  Q*°( J  +  *m  -  rm>  r) 
r=0 

Since  the  curves  /0,...,/p  are  distinct,  the  a's  corresponding  to  different  curves  are  different 
from  each  other.  Therefore  the  coefficients 

{a(/i,/2)|/i  +  ml 2  =  mi  +  j;  } 
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can  be  uniquely  found  by  solving  the  above  Van  der  Monde  system  of  equations.  The  same 
procedure  can  be  applied  for  finding 

{a(*i>*2)l*i  +  mh  =  n(m  +  1)  -  im  -  j\  } 


from  the  quantities 

Therefore  we  have  shown  that  in  step  i  +  1  we  can  uniquely  determine  the  quantities 

W*i)  *2)|*i  +  »™*2  =  rm, m(i  +  1)  -  1,  n(m  +  1)  -  («  +  l)m  +  1, n(m  +  1)  -  im;  } 

This  completes  the  induction  and  proves  the  theorem.  It  is  worthwhile  to  mention  that  since 
in  the  fth  step  of  the  algorithm  we  find  2mi  coefficients  of  p(w,  2),  for  arbitrary  values  of  n  and 
m,  there  is  no  guarantee  that  the  number  of  interpolation  points  is  equal  to  the  number  of  the 
unknown  coefficients.  In  fact  unless  there  is  an  integer  p  such  that 

p 

y^(m  +  l)n  -  2mi  +  1  =  (n  +  l)2 
1=0 


we  need  more  data  points  than  unknown  coefficients. 


Appendix  B 

A  Bound  on  the  Number  of  Finite 
Common  Zeros  Based  on 
Polynomial  Degree 


Recall  that  Bezout’s  theorem  is  concerned  with  determining  the  number  of  common 
zeros  of  two  bivariate  polynomials.  It  is  restated  here  for  convenience, 

Theorem  B.l  If  p(w,z)  and  q(w,z)  are  bivariate  polynomials  of  total  degree  r  and  s  respec¬ 
tively  with  no  common  factors,  then  there  are  at  most  rs  distinct  pairs  (w,z)  where 

p(w,z)  =  q(vu,z)  =  0  (B.l) 

Since  Bezout’s  theorem  is  concerned  with  total  degree  as  opposed  to  degree  in  each  variable, 
it  pertains  most  generally  to  polynomials  whose  coefficients  have  triangular  support  as  shown 
in  Figure  B.l. 

In  our  case  of  reconstruction  from  level  crossings,  this  corresponds  to  an  image  which  has  a 
triangular  support.  On  the  other  hand,  many  times  one  is  interested  in  images  with  square  or 
rectangular  support  as  shown  in  Figure  B.2. 


For  the  case  when  the  polynomials  under  consideration  have  rectangular  support,  we  are 
able  to  lower  the  bound  on  the  number  of  common  finite  zeros  from  the  bound  set  by  Bezout’s 
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Figure  B.l:  Non-zero  coefficients  for  a  polynomial  of  total  degree  2. 
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Figure  B.2:  Non-zero  coefficients  for  a  polynomial  of  degree  (2,2). 
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theorem.  Specifically  if  the  two  relatively  prime  polynomials  p(w,  z)  and  q(w,  z)  are  given  by 


P(«>, 


M„  N„ 

*)  =  £  E 

m=0  n=Q 


(B.2) 


and 


M,  N< 

£  9m,n^m2n 

m= 0  n=0 


(B.3) 


the  upper  bound  on  the  number  of  common  finite  zeros  set  by  Bezout’s  theorem  is  (Np  + 
Mp)(Nq  +  Mq).  Our  objective  is  to  establish  a  tighter  upper  bound  on  the  number  of  common 
finite  zeros  of  p(ti/,z)  and  ^(ti/jz).1 


Before  proceeding,  we  need  to  review  several  results  concerning  the  resultant  of  polynomials 
in  one  or  two  variables.  The  resultant  Rpq  of  two  one-dimensional  polynomials  p(w)  and  q(w) 


n=0 


(B.6) 


is  defined  [37]  as  the  determinant  of  the  (Mp  +  Mq)  by  (Mp  +  M?)  matrix 


Po  Pi 
0  Po  Pi 


PM„  0  0  .  .  0 

PM,,- 1  PMf  0  .  .  0 


0  0 
7o  9i 
0  qo  ?i 


7M,  0  0 

9M,-i  7M,  0 


0 


PMp 

0 

0 


(B.7) 
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‘The  derivation  presented  here 
described  in  [39]. 
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is  due  to  a  collaboration  between  D.  Izraelivitz  and  the  author  and  is  also 


A  basic  property  of  resultants  is  stated  in  the  following  theorem  [37], 


Theorem  B.2  When  the  polynomials  p(ui)  and  q(w)  have  numerical  coefficients,  a  necessary 
and  sufficient  condition  that  they  shall  have  a  finite  or  infinite  common  root  is  that  Rpq  =  0. 

Consider  now  the  two  relatively  prime  bivariate  polynomials  p(w,  z)  and  q(w ,  z),  defined  in 
(B.2)  and  (B.3),  expressed  as  polynomials  in  w,  with  coefficients  which  are  polynomials  in  z, 

p{w,  z)  =  po(z)  +  pi(z)w  H - +  pmv{z)wMv  (B.8) 

g(u>,  z)  =  q0(z)  +  qi[z)w  +  •  •  •  +  qMii{z)wM<>  (B.9) 

We  can  define  the  resultant  of  p(w,  z)  and  q{w,  z)  with  respect  to  u;  as  the  determinant  of  the 
(Mp  +  Mq)  by  (Mp  +  Mq)  matrix,  M(z),  with  polynomial  entries: 


M(z)  = 


'  Po(z) 

Pi(*) 

PMr{z)  0 

0 

0 

0 

Po(z) 

Pi(z)  ■ 

PM,-l(z)  PMp(z) 

0 

0 

0 

0 

■  Po(z ) 

■  PMv{z) 

90(2) 

9i(*) 

9m,  (2) 

0 

0 

0 

0 

9o(z) 

9i(^)  • 

9M,-l(2) 

9M„(2) 

0 

0 

0 

0 

0 

90(2) 

•  9M,  (2)  . 

(B  10) 

This  resultant  is  a  function  of  the  remaining  variable  z  and  is  denoted  by  Rpq(z).  Expanding 
the  determinant  of  the  above  matrix  and  taking  into  account  that  each  Pi{z)  and  qt{z)  is  of 
degree  at  most  Np  and  respectively,  we  can  conclude  that  Rpq(z)  is  a  polynomial  of  degree 
NpMq  +  NqMp  or  less.  It  can  be  shown  [38],  that  if  p(u>,  z)  and  q(w ,  s)  are  relatively  prime  then 
Rpq(z)  is  not  identically  zero.  Moreover,  if  (ui0,  zq)  is  a  common  zero  of  p(w,  z)  and  q(w,  z)  then 


Rpq(zo)  =  0-  Thus,  the  zero  sets  of  p(w,  z)  and  q( w,  z)  have  at  most  NpMq  +  NqMp  values  of  z 
in  common. 

As  our  argument  stands,  we  have  not  yet  placed  any  tight  limit  on  the  number  of  intersec¬ 
tions  of  p(w,z)  and  q(w,z)  since  for  each  z,  which  is  a  root  of  Rpq{z)  there  could  be  a  large 
number  of  u>3  ,  such  that  for  each  j, 

p(wj,Zi)  -  q{wj,Zi)  -  0  (B  11) 

In  order  to  specify  the  number  of  uij  for  each  z<,  we  need  to  study  the  behavior  of  Rpqiz) 
in  the  vicinity  of  each  zt . 

Theorem  B.3  If  at  each,  zo  there  are  k  values  of  w,  uij,  such  that 

p{wj,z0)  =  q(wj,z0)  —  0  for  j  —  1,  •  •  • ,  A:  (B.12) 

then  Rpq(z)  has  a  zero  of  multiplicity  at  least  k  at  zq. 

The  above  theorem  implies  that  p(w,z)  and  q[w,z)  as  defined  in  equation  (B.2)  and  (B.3) 
have  at  most  NpMq  -I-  NqMp  zeros  in  common. 

In  order  to  prove  Theorem  B.3  we  need  to  review  some  results  on  matrices  with  polynomial 
entries,  relating  to  the  Smith  Normal  Form  [56], 

Theorem  B.4  Let  A(x)  be  an  n  by  n  polynomial  matrix  of  normal  rank  r.  We  can  find 
unimodular  matrices  {P{x),Q{x)} ,  such  that 

B(x)  =  P(x)A<x)Q(x)  (B.13) 

and 

1.  B(x)  is  a  diagonal  polynomial  matrix  called  the  Smith  Normal  Form  of  A(x). 

2.  The  first  r  diagonal  elements  of  B(x)  are  monic  polynomials  pi(x),  p i(x),  ■■■,  Pr{z). 

3.  The  remaining  diagonal  elements,  if  any,  are  zero. 
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f.  Pi(x)  divides  pi+i(x)  fori  =  l,---,r-  1 


Normal  rank  in  the  above  theorem  is  defined  in  [56] .  The  unimoduiar  polynomial  matrices 
of  the  above  theorem  are  defined  to  have  nonzero  constant  determinant  independent  of  x. 
Therefore,  if  r  =  n,  i.e.,  if  A(x)  has  full  normal  rank  then, 

det(£(x))  =  p;(x)  (B-14) 

«=i 

Also,  from  part  (4)  of  Theorem  B.4  we  can  conclude  that  if  p,(x)  =  0  then  p*( x)  =  0  for 
k  >  i.  From  the  above  theorem,  we  can  derive  the  following, 

Theorem  B.5  Let  A(x)  be  a  polynomial  matrix  of  full  normal  rank.  If  A(x o)  has  rank  defi¬ 
ciency  of  k  then  det(A(z)),  has  a  zero  of  multiplicity  at  least  k  at  x o- 


Proof: 

Using  Theorem  B.4  we  can  find  B(z),  the  Smith  normal  form  of  A(x).  Since  P(x)  and  Q{x) 
are  unimoduiar,  B(x)  is  of  full  normal  rank.  Furthermore,  the  ranks  of  B(x)  and  A(x)  at  each 
value  of  x,  including  xo,  are  equal.  Therefore  B(xo )  has  rank  deficiency  of  k.  This  means  that 
pn(x0)  =  pn_l(z0)  =  •••  =  p„_fc+i(z0)  =  0.  Therefore,  from  (B.14)  det(B(x))  has  a  zero  of 
multiplicity  at  least  k  at  zo-  Since  the  determinant  of  A(x)  is  within  a  constant  factor  of  that 
of  B(x),  A{x)  also  has  a  zero  of  multiplicity  at  least  k  at  xq. 

□ 


Using  Theorem  B.5,  we  can  return  to  Theorem  B.3, 


Proof : 


(Theorem  B.3)  Suppose  that  for  zq  there  are  k  common  finite  zeros  Wj,  j  —  1,  •  •  • ,  k  between 
p(u/,  z)  and  q(w,  z).  Then  the  matrix  M(zo)  defined  by  (B.10)  must  have  k  linearly  independent 
null  vectors  given  by: 

|1  (B.  15) 


for  j  =  1, ■••,k.  Thus  M(zq)  has  rank  deficiency  of  k.  Furthermore,  since  p(w,z)  and  q(w,z) 
have  no  common  factors,  fipq(z)  is  not  identically  zero,  so  M(z)  has  full  normal  rank.  From 
Theorem  B.5  then,  Rpq(z)  has  a  zero  of  multiplicity  of  at  least  k  at  zo 


%'t  >% ; t »%/ 


From  Theorem  B.3,  and  the  fact  that  Rpq{z)  is  a  polynomial  of  degree  NpMq  +  NqMp,  we 
get  immediately,  the  main  result  of  this  appendix, 


Theorem  B.6  Let  p(w,z )  and  q(w,z)  be  two  polynomials  of  degree  (Mp,  Np)  and  ( Mq,Nq ) 
with  no  common  factors,  then  there  are  at  most  NpMq  +  NqMp  pairs  of  finite  complex  numbers 
(w,z)  such  that 

p(w,z)  =  q(w,z)  =  0  (B-16) 


Appendix  C 

A  Result  in  Algebraic  Geometry 


In  this  appendix,  we  will  briefly  go  over  few  definitions,  and  then  prove  a  result  in  algebraic 
geometry,  which  is  used  in  the  proof  of  Theorem  (4.2). 

Consider  the  polynomials 


pt(x0,...,xn)  =  0  ,  t  =  1, ..., r 


defining  an  irreducible  algebraic  set  V  over  K.  Let  us  denote  the  r  x  n  Jacobian  matrix 
associated  with  V  by  (|^*-).  Points  x t  of  V  at  which  the  Jacobian  matrix  assumes  is  maximum 
rank  are  called  ordinary  points.  Any  point  of  V  which  is  not  ordinary  is  said  to  be  singular 
(57j.  If  the  singular  points  of  V  are  removed,  we  obtain  a  manifold  whose  dimension  defines 
the  topological  dimension  of  V  *.  The  dimension  of  a  reducible  algebraic  set  is  defined  to  be 
the  maximum  of  topological  dimensions  of  its  various  irreducible  components. 

The  complex  topological  dimension  of  V  C  Cn,  which  is  defined  to  be  the  dimension  of  its 
associated  complex  manifold,  is  also  given  by  [59|: 


,  dU  , 


‘Precise  definition  of  manifolds  is  included  in  [ 58 1 


where  P  ranges  over  the  points  in  V.  By  definition,  the  real  dimension  of  V  is  twice  its  complex 


dimension.  For  convenience,  we  will  use  the  term  topological  dimension  to  denote  real  dimension 
unless  specificed  otherwise.  If  V(R)  denotes  the  real  part  of  V  c  Cn,  we  have  [60]: 

dimtopV(R)  <  i dimtopV 

where  dimtop  denotes  the  topological  dimension.  V(R)  is  said  to  be  of  maximal  topological 
dimension  if  its  dimension  is  exactly  half  the  dimension  of  V. 

We  are  now  ready  to  prove  the  following  theorem  2  which  is  ultimately  used  in  the  proof  of 
Theorem  (4.2). 


Theorem  C.l  Let  V  c  CN  be  the  set  of  complex  zeros  of  polynomials  Then  if  V 

is  irreducible  and  V(R),  the  real  points  ofV  has  maximal  topological  dimension,  then  V (R)  is 
Zariski  dense  in  V.  That  is,  every  polynomial  that  vanishes  on  V  (R)  must  vanish  on  all  ofV. 


Proof: 

We  will  prove  the  above  result  by  contradiction.  If  a  polynomial  /  vanishes  on  V(jR)  and 
does  not  vanish  on  V ,  then  let  W  C  V  be  the  set  of  real  and  complex  zeros  of  /  which  are  in 
V .  Since  V  is  irreducible,  and  W  is  a  proper  subset  of  V, 

dimtopW  <  dimtopV  (C.l) 


Since  the  real  part  of  W  and  V  are  the  same,  then 


dimtopW{R)  =  dimtopV(R) 


Using  the  above  equation  and  the  assumption  of  maximal  topological  dimension  for  V{R),  we 
get 


dimtopW(R) 

Furthermore,  as  we  mentioned  earlier 

dimtopW  (R) 


^  dimtopV 

(C.2) 

dimtopW 

(C.3) 

2This  theorem  and  it*  proof  were  suggested  by  Prof.  M.  Artin  at  MIT. 


From  equations  (C.2)  and  (C.3)  we  conclude  that: 


dimtopV  <  dimtopW 


which  contradicts  inequality  (C.l). 


VClVJv*  V_ s'  v.'V \-V 
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